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Abstract 

The terms renormalization and renormalization group arc explained by reference 
to various physical systems. The extension of renormalization group to turbulence 
is then discussed; first as a comprehensive review and second concentrating on the 
technical details of a few selected approaches. We conclude with a discussion of the 
relevance and application of renormalization group to turbulence modelling. 
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Part I. Introduction 


Part I 

Introduction 

1 Preliminary 

The term “RG” (or “RNG”) has now become familiar to the turbulence modeling 
and computational fluid dynamics (CFD) community. However, understanding just 
what RG is actually about is often viewed as an obscure mathematical problem. 
Accordingly, we feel that there is a need for an appropriate article for this diverse 
community. The present review article, which surveys many aspects of applications 
of RG to turbulence, is intended to fill this gap. 

Turbulence in fluids has presented a notoriously difficult problem for more than a 
century. As discussed in standard text books on fluid motion, the major difficulty is 
the simultaneous existence of many different space and time scales at high Reynolds 
number. All these scales are of equal importance, hence the theoretical difficulty. 

Now, RG was a tool originally developed for dealing with such problems. Indeed, 
the work of Wilson in developing the RG method for critical phenomena led to his 
being awarded the Nobel prize in Physics in 1982. In The Nobel Prize Winners: 
Physics (Magill 1989), we read: 

“Wilson clarified a natural phenomenon that has puzzled scientists throughout the 
centuries: the behavior of substances at critical points, or phase transitions. With his 
renormalization group theory, Wilson divided this seemingly insoluble problem into 
a number of small ones.” 

1.1 General remarks 


As we have said, at high Reynolds numbers ( R e ), a wide range of length and time 
scales in turbulent flows become excited. Direct numerical simulations (DNS) which 
are three-dimensional, time-dependent numerical solutions in which all significant 
scales of motion are computed without any modeling are currently restricted to low 
turbulence Reynolds numbers and simple geometries (Reynolds 1990). Some recent 
DNS of homogeneous and isotropic turbulence at 512 3 resolutions are Chen et al. 
(1993), Wang et al. (1996), and Yeung and Zhou (1997). 

Although DNS can be extremely useful in many areas related to the study of 
turbulence physics and the assessment of theories (ibid), virtually all scientific and 
engineering calculations of nontrivial turbulent flows, at high R e , are based on some 
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type of modeling (Speziale 1991). Most modeling efforts have been focused on subgrid 
scale models for large-eddy simulations and Reynolds stress models the two forms 
of modeling that will be discussed in this article. 

Large-eddy simulations (LES) are a logical ‘modeling extension’ of DNS. They 
arc based on the observation that the small scales are more universal in character 
than the large, energy containing scales of motion. In LES, the three-dimensional 
time-dependent motion of these large scales are computed directly while the effects of 
the small scales on the large scales are modeled. This procedure leads to the so called 
‘subgrid scale modeling’ problem; see Rogallo and Moin (1984), (Reynolds 1990), 
Lesicur and Mctail (1996), and Moin (1997). The traditional concept of eddy viscosity 
(Smagorinsky, 1963) is that the subgrid turbulence simply enhances the viscous term 
in the momentum equation, and that the value of the enhanced eddy viscosity can 
be determined by equating the resulting increase in the energy dissipation rate to the 
rate at which energy is transferred from resolved to subgrid scales. LES is useful in 
the study of turbulence physics at high Reynolds numbers that are unattainable by 
DNS. 

LES is also intended to be useful in the development of turbulence models for 
the prediction of the complex flows of technical interest where simpler modeling ap- 
proaches fail (Reynolds 1990). However, all applications of LES have been in simple 
geometries where Van Driest damping could be used an empirical approach that 
generally does not work well when there is flow separation (Speziale 1991). 

We note that LES and DNS are recent developments, being direct products of 
the rapid increase in computer capabilities. The Reynolds stress modeling concept, 
however, was introduced by Osborne Reynolds as long ago as 1895. Reynolds stress 
models remain the most important tool for the highly repetitive engineering calcula- 
tions associated with design (Reynolds 1990). The goal of Reynolds stress modeling 
is to provide reliable information about first and second one-point moments (e.g., the 
mean velocity, mean pressure, and turbulence intensity), which is usually all that is 
needed for design purposes (Speziale 1991). In principle, this approach attempts to 
model all of the fluctuating scales. 

Recently, the methodology of RG has attracted considerable attention as a sys- 
tematic approach to subgrid scale modeling in LES and Reynolds stress modeling. In 
principle, the RG removal of only the smallest scales generates a subgrid scale model 
while Reynolds stress models are obtained in the limit when all fluctuating scales are 
removed. This point of view is due to Yakhot and Orszag (1986). 

1.2 Overview 

Because we are addressing such a diverse community, some innovation in the structure 
of the review is indicated. Part II is intended for those readers who wish to obtain 
a basic knowledge of “RG” and what are the main results from this approach to 
studying turbulence; and in that sense it is complete within itself. 
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Part III is for those readers who would like to go more deeply into RG theory. 
Here we present the detailed mathematical procedures of two RG approaches. In 
particular, we clarify their major approximations for the benefit of the reader who is 
willing to pursue this line of research. In this part, we have included several related 
pieces of work and have also attempted to compile a comprehensive reference list. 

Some turbulence modelers or CFD workers may wish to skip Part III and proceed 
directly to Part IV. Here we present recent efforts to develop turbulent Reynolds 
stress models based on RG. 
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Part II 

General description of RG 

2 What Is Renormalization Group? 

In this section we shall explain something of the origins and meaning of RG. In the 
process we shall define terms and introduce concepts which will be needed when we 
consider the question of how RG can be applied to fluid turbulence. 

We should begin by noting the difference between the concept of renormalization 
and renormalization group. The renormalization procedure was introduced originally 
as a method of removing divergences in quantum field theory. The renormaliza- 
tion group theory, on the other hand, was used for improving perturbation theory 
(Gell-Mann and Low, 1953) by exploiting the non-uniqueness in the renormalization 
procedure (Stuechelberg and Petermann, 1953). We can illustrate this idea by con- 
sidering a problem which may be discussed entirely in terms of classical physics. This 
is the electron gas, which consists of a gas of electrons in a uniform background of 
positive charge, so that the system is, overall, electrically neutral. Such a gas can be 
a model for the conduction electrons in a metal, for a plasma, or for an electrolyte. 
This allows us to introduce the general idea of renormalization, and is the subject of 
Section 2.1. In order to introduce the more specific topic of renormalization group, 
we consider as an example the phase transition from para- to ferromagnetism. This is 
dealt with in Section 2.2. We will review the applications of RG to critical phenomena 
in Sections 2.3. 

2.1 Renormalization 

It is a cardinal assumption of physics that the macroscopic properties of a material, 
such as its dielectric constant or elastic modulus, will not in general depend on the 
size of the sample of material under consideration. So, for a microscopic theory, t his 
implies that we must be able to take limits, corresponding to the distance between 
constituent microscopic particles shrinking to zero, or to the overall size of the system 
going to infinity, without upsetting our theory. However, if we take our example of 
an electron gas, it is easily seen that problems may arise in both these limits. 

Let us suppose, for instance, that we wished to work out the total electrostatic 
potential energy of all the electrons in the system and that we did that naively by 
adding up the contributions from individual electrons, each electron contributing the 
usual Coulomb ‘1/r’. It is obvious that a continuum limit would require r — > 0, with 
each potential diverging in that limit. Less obviously, the limit of infinite system size 
also leads to a divergence. Suppose we take the test charge to be at the origin and 
assume that the surrounding charges have a uniform density p, say. Then, as the 
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number of electrons in a spherical shell between r and r + dr goes as pr 2 dr, it is easily 
seen that this overpowers the Coulomb potential to give a total potential which varies 
as 2n peR 2 , where e is the electronic charge and R is the system size. 

The first of these divergences is usually referred to as being ‘ultra-violet’, because 
short distances correspond to high spatial frequencies. Similarly, the long-range di- 
vergence in the second case, is referred to as being an ‘infra-red divergence’. From our 
point of view, the second case is the more interesting. In a charge neutral plasma, the 
energy at a point does not diverge with system size because a cloud of ions around 
each individual electron acts to screen the long-range Coulomb potential and turns 
it into an exponentially decaying potential (for stationary particles). Theoretically, 
this effect shows up when we take the interactions between individual electrons into 
account. In perturbation theory it would correspond to summing a certain class of 
terms to all orders. In fact, the potential at lowest nontrivial order, the potential 
becomes ecxp[— r/lo\/r, where Id is the Debye length and depends on the density of 
electrons. 

This is known as a screened potential. But, an alternative interpretation is to say 
that we replace the electron charge e by the renormalized electron charge e exp[— r/lp}- 
This replacement of a ‘bare’ particle by a ‘dressed’ particle, in order to take account 
of the effect of interactions, is an important technique in many-body physics and 
nowadays is widely referred to as renormalization. It should be noted that a charac- 
teristic feature of the process is that the renormalized properties (as in our example) 
are scale-dependent. 

This discussion serves to introduce the idea of renormalization. Our next step is 
to consider what is meant by renormalization group , and it turns out that this is just 
one way of achieving renormalization. 


2.2 Renormalization group 

RG was based on the obvious equivalence between taking the continuum limit and 
making a scaling transformation of the basic variables such as position or momentum. 
However, this theory was regarded as rather opaque in terms of its underlying physics, 
and from our point of view, it is better to consider its more modern incarnation in the 
theory of critical phenomena, in which it is seen as a way of dealing with problems 
involving a large range of length or time scales. Its possible relevance to turbulence 
might now begin to emerge; but we shall first consider another simple physical system: 
the lattice-spin model of a ferromagnet. 

Magnetism arises because spins at different lattice sites tend to become aligned. 
This tendency is opposed by thermal effects, which promote random orientation of 
the spins. At any fixed temperature, alignments occur on length scales ranging from 
the lattice spacing Lq up to some correlation length £. Take a temperature T 0 which 
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is sufficiently high that £(T 0 ) is of order L 0 . It is evident that £ will increase as T 
decreases below T 0 . In fact, £ actually becomes infinite as T approaches the Curie 
temperature , which defines the critical point for this problem. At this critical point, 
fluctuations on all scales from Lq up to the size of the material specimen can occur, 
and therefore a net overall magnetization can appear. 

The theoretical problem is to calculate the partition function and thence the 
thermodynamic properties of the material. The partition function is a sum over all 
possible spin configurations; although this sum is easily evaluated if the spins at 
different sites are independent, coupling between spins makes this problem difficult. 
The fact that all length-scales are (in principle) equally important makes it difficult 
to know what to do. 

The problem is eased to some extent by considering models, such as the Ising 
model, in which spin vectors are Boolean in character and are either ‘up’ or ‘down’, 
with no intermediate states permitted. Another helpful feature of this Ising model 
is that the spins only have nearest-neighbour interactions. Then the application of 
RG to this problem can be understood as a form of coarse-graining, in the following 
sense. The version of the RG method, called the block-spin transformation, consists 
in breaking down an intractable problem with multiple scales of length into a sequence 
of smaller problems, each of which is confined to a single length scale. This procedure 
consists of three steps. First, the lattice is divided into blocks of a few spins each 
(N 2 , say, for a two-dimensional lattice). Then each block is replaced by a single spin 
whose value is the average of all the spins in the block; there the average is determined 
by the majority rule. Therefore, a new lattice is created with N times the original 
spacing and 1/A the density of the spins. Finally, the original scale is restored by 
reducing all dimensions by a factor of N. 

We now provide a brief description on the RG transformation with the block 
spin method. Our starting point is the Hamiltonian H 0 , associated with two spins 
separated by a distance L 0 (i.c. the lattice spacing). Then we calculate an effective 
Hamiltonian H\, associated with regions of size 2 L 0 , by averaging out the effects of 
scales L 0 . Next we calculate H 2 , associated with regions of size 4L 0 , with the effects 
of scales less than or equal to 2 L 0 averaged out. This well known block spin method 
is due to Kadanoff (1977). 

The above process can be expressed in terms of a transformation 7~, which is 
applied iteratively: TH 0 — > ► H\,TH\ — > H 2 , TH 2 — ► . . .. At each stage, the length 

scales are changed: L 0 — ► 2 L 0 , 2Lq — > 4Lo . . ., and in order to compensate, the spin 
variables are also scaled in an appropriate fashion (sec next subsection). It is this 
rescaling which leads to renormalization, and the set of transformations {T} defines 
a semi- group: hence renormalization group. If iterating this transformation leads to 
the result that an integer N exists such that 

T H n = H n+X = H n , for alln > N 

then H n = H N (say) is referred to as a fixed point. In the case of critical phenomena, 
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this corresponds to a critical point. We conclude this brief discussion of the magnetic 
case by indicating how its methodology may be made to come closer to a continuum 
problem (as in fluid motion), which is what concerns us. 

Suppose we make the Fourier transformation from x-spacc to wavevector k-space. 
Further suppose that we have eliminated (in some way) the smallest length scales, less 
than or equal to A -1 , say, so that we arc left with a Hamiltonian which is a function 
of continuous variables in k-space, H{ A). Then we may formally eliminate further 
modes in a band A/6 < A; < A, as follows (1 < b < oo). Denote the three-dimensional 
spin held by ^(k), where a = 1, 2, or 3, and take the following steps: 

1. Integrate out all .S Q (k) for which A/6 < k < A. 

2. Rescale the remaining modes of the spin held by enlarging wavevcctors by a 
factor 6. 

3. Multiply each S Q (k) by a constant factor Q. 

The parameter 6 is known as the spatial rescaling factor, while Q, is the spin rescaling 
factor, and involves technicalities of the magnetic case, which we shall not pursue 
here. If we write the above three-step process as 

H' = T b H, 

then the set %, 1 < 6 < oo is called the renormalization group. It is often remarked 
that, since T b 1 does not exist, this is strictly a semi-group. 

We shall now consider RG as applied to the 2D Ising model. 

2.3 Renormalization group theory in critical phenomena 

2.3.1 Iterative Hamiltonian structure 

As we discussed already, an important and early application of RG was to critical 
phenomena. Here we shall provide a more detailed review of the calculation of Wilson 
(1975; Wilson and Kogut, 1974) for the 2D Ising model. In the 2D Ising model 
of ferromagnetism, one considers a lattice of points and on these lattice sites are 
particles (called ’spin’ particles) which (i) can take on only two values [ i.e., ‘spin up’ 
or ‘spin down’], and (ii) can only interact with its nearest neighbors. The system is in 
thermal equilibrium, with a lattice temperature kjfT . For simplicity, wc shall consider 
a square lattice with many spin particles. ^From equilibrium statistical mechanics 
(Huang, 1963), it is known that all thermodynamic properties of this system can be 
determined from the partition function 

Z = J2 ex P H o[s] (1) 
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where the Hamiltonian H a , for the 2D Ising model, 

H 0 [s] = K S nSn+i (2) 

K = — J/k B T, where J is the spin-spin interaction coupling constant [with J > 0 for 
ferromagnetism]. The sum in (1) is evaluated over all possible spin configurations. 
The outer sum in (2) is over all lattice points, while the inner sum, for the represen- 
tative lattice site vector n = (n i, r i 2 ) ■ is a sum over the four nearest neighbors of the 
square lattice 

s n s n+i — s n s (m + l,ri 2 ) T S n S( ni _x n2 ) + S n S( nj n2+ j) + S n S( ni (3) 


The difficuly in determining the thermodynamic properties of this microscopic 
Ising model is in the evaluation of the summation in (1) : there are just too many 
degrees of freedom 2 . The important things to note are (i) the original problem in- 
volved only nearest-neighbor interactions, (ii) the coupling strength K is related to 
the strength of the interactions, J, and the lattice temperature /c fl T. 

2.3.2 Iterative decimation of the number of degrees of freedom 

Suppose that the degrees of freedom are decimated by restricting the summation in 
Eq. (1) to run only over half the lattice sites. Thus the partition function Z must 
now be determined from 

Z = Y.e*p[HM\ ( 4 ) 

1 

where the new lattice sites are denoted as er, the sum is over these sites only, and H\ 
is a Hamiltonian to be determined so that (4) yields the identical partition function 
as Eq. (1). 

It can be shown that the new Hamiltonian H\ [a] takes the form 


Hi [a] = (N 2 /2)A + 2BJ2 ViVk + + C ]T] (5) 

nn nnn p 

where “nn” denotes nearest-neighbor coupling, “nnn” next-nearest-neighbor coupling 
and l p’ the plaquette four-spin interactions (Hu, 1982) with 

A(k) = ln2 + (1 / 8)[ln(cosh4K) 4- Aln(cosh2K )] (6) 

B(k) = (. \/%)ln{cosMK ) (7) 

C{k) = ( 1/8) [ln(cosMK) - Aln(cosh2K)\ (8) 

2 Wc are actually using this 2D Ising model as an example of how to treat problems that have 


so many degrees of freedom that new techniques must be introduced to solve them. The 2D Ising 
model is actually exactly soluble, but that is not of direct concern to us. 
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Here A is simply a normalization. Eq. (5) should be compared to the original Ising 
Hamiltonian H 0 , eq. (2). First we see the expected changes: (a) the nearest neighbor 
interactions term has strength 2B(K), rather than K, and (b) there is a normaliza- 
tion constant A(K) introduced due to the change in the correlation length of the 
decimated system. What is of significant interest is the appearance of new, non- 
Ising, interactions: diagonal nearest neighbor interactions [strength B(k)], and four 
spin interactions [with strength C(k)]. Therefore, although the original Hamiltonian 
has only nearest-neighbor coupling, other types of interactions arc generated by the 
renormalization group to compensate for the reduction in the number of degrees of 
freedom. This is a general feature of the renormalization group (Hu, 1982). 

We now proceed to the second iteration to determine the Hamiltonian ff 2 for 
N 2 / 2 2 lattice sites. H 2 is to be determined so that the partition function Z remains 
invariant under decimation of spins. If all the 4 terms are kept in (5) in determining 
H 2 , then one cannot calculate the sum over the decimated spins exactly. To pro- 
ceed, one can start with the nearest-neighbor interaction and treat the other types 
of interactions perturbatively. Wilson proceeds under the assumption that the 4-spin 
interaction (and any higher order spin interaction) is weak and can be neglected, but 
it is critical that the U uuu n interactions be retained. The self-consistency of these 
assumptions can be checked at the end of the calculation. 

Hence, if at the i - I th stage, the Hamiltonian is H t _i[s], with its dependence 
on the Ising (original) nearest-neighbor interaction strength K^\ and the (non-Ising) 
diagonal-ncarcst-neighbor interaction L*_i, then at the i th iteration, the Hamiltonian 
Hi[ s] is a function of K x and L x . It is important to note that only in degenerate cases 
is the Hamiltonian form invariant under RG transformations. 

2.3.3 RG transformations 

Now for sufficiently weak interactions, one can expand the interaction strength coef- 
ficients in a Taylor scries in K — — J/KbT : 

B{K) = ln{coshAK)/8 « K 2 + 0(K 4 ) (9) 

C(K) = ( 1 / 8)\ln{coshAK ) — Aln(cosh2K )] = 0(K 4 ), (10) 

allowing us to drop the 4-spin interaction term in (5). Note further that the coeffi- 
cients for “u«” and “ituu” interactions are of the same order in K, requiring us to 
keep these newly created diagonal “«uw” interactions. To determine the RG trans- 
formation, we must note the following: (a) the diagonal nearest-neighbor interaction 
at the i th stage, L t , arises from the coupling coefficient = Kf_ 1 + O ( K 4 _ x ) , 

on using (7): i.c., 

Li = K?_ i (11) 

(b) the nearest- neighbor interaction at the i th stage arises from 2 terms: (i) the 

nearest-neighbor interactions at the i — 1 th stage, with coupling coefficient Z?(/\,_i) = 
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K?_ j , and (ii) the diagonal-nearest-neighbor interactions at the i — 1 th stage which 


yields a coupling coefficient L t _i to K x . This is because the diagonal- nearest-neighbor 
interaction at the i — 1 th stage will become the nearest- neighbor interaction at the 
i th stage due to the decimation of spins. 

Hence the RG transformations to proceed from the i — l t/l to the i th stage are 

K x = 2 KU + U - 1 

(12) 

U = Kl 1 

(13) 

with 

L 0 = 0 

(14) 

since initially we are dealing with an Ising model which, by definition only has nearest- 
neighbor interactions. 

The fixed points, AT* and L*, of (12)- (14) are immediately determined from: 

K< = 2 Kl + L, 

(15) 

U = Kl 

(16) 

There arc 3 possible fixed points: 


K + 1 ~ 0, L*i = 0 

(17) 

K* 2 = OO, L+2 — OC. 

(18) 

K * 3 - 1/3, u 3 = 1/9 

(19) 


and the initial domain of attraction in the K-L plane for these fixed points are readily 
determined (with L 0 = 0): (a) for 0 < K 0 < 0.39209, L 0 ~ 0 : stable fixed point 
is (0,0); (b) for K 0 > 0.3921, L 0 = 0 : stable fixed point is (oo,oo); (c) for K 0 = 
0.39209, Lq = 0 : unstable fixed point is (1/3, 1/9) 

In the Ising problem, the relevant fixed point is the saddle point (c). The fixed 
point (a) corresponds to K * = 0 = J/ksT, » , i.e., T* — ► oo. Physically, this corresponds 
to a perfectly disordered system, with only short scale fluctuations that tend to zero 
as the number of RG iterations. The fixed point (b) corresponds to the temperature 
T, = 0, and this corresponds to a perfectly ordered system, again with negligible 
fluctuations. Hence, as far as critical phenomena are of interest, these fixed points 
are unimportant. However, the fixed point (c) is important since even after all the 
RG transformations, there arc still many length scales present. Critical exponents 
prove to be properties of this fixed point. 
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3 Application of RG To Turbulence 


3.1 Similarity and differences between fluids and critical phe- 
nomena 

The macroscopic motion of a fluid may be thought of as exhibiting two ‘phase tran- 
sitions’. That is, as wc increase the Reynolds number of a given flow, we may expect 
first a transition to turbulence and second a transition to scaling behaviour of fully 
developed turbulence. In the latter case, this means the existence of an intermediate 
range of wavenumbers in which the energy spectrum takes the form of a power-law. 
It is this transition which concerns us here. Specifically, we will assume the iner- 
tial range ‘5/3’ law for the energy spectrum developed by Kolmogorov [1941]. As 
an example, consider a stationary homogeneous isotropic turbulence where energy 
is injected with a forcing spectrum peaked a wavenumber k 0 (= 1 /Iq), where l 0 is 
the integral length scale. The energetics of eddies with wavenumbers smaller than 
a given wavenumber k are determined by an equilibrium between the injection, the 
losses due to viscous dissipation and the energy flux function II (k) that moves energy 
to higher wavenumbers (i.c., to smaller scales). Here, 1 1 (A:) depends upon the inertial 
terms in the Navier-Stokes equation, and is expressible in terms of the triple veloc- 
ity correlation (see for example, Rose and Sulcm, 1978). The problem is defined by 
three parameters plus the underlying Navier-Stokes equation. These parameters arc 
the kinematic viscosity, the total rate of dissipation, e, which equals the energy flux 
function, and the characteristic lentil / 0 . To make contact with the critical phenom- 
ena, the external length scale Iq might be considered analogous to a lattice constant 
(Nelkin, 1974). 

It is of course difficult to identify the onset of scaling behaviour as such and a con- 
venient way to think about this problem is to consider the case where the Reynolds 
number (based on the turbulent microscale, say) is large enough for a power-law spec- 
trum to exist over an appreciable range of wavenumbers. Now define a local Reynolds 
number, using the inverse of wavenumber as the length scale. Then as we scan through 
the wavenumbers, from the maximum towards smaller values, wc arc, in effect, in- 
creasing the local Reynolds number to the point where power-law behaviour begins 
(In other words, the boundary between the inertial and viscous ranges). Conceptu- 
ally, therefore, the simplest application of RG to turbulence involves the progressive 
elimination of high- A; modes, with the fixed point corresponding to a renormalized 
viscosity in the inertial range (The equivalent, in terms of CFD, would be a subgrid 
effective viscosity, in large-eddy simulations). Wc shall carry out these operations on 
the Navier-Stokes equations. 
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3.2 RG applied to the Navier-Stokes equations 

Wc introduce the Fourier components of the velocity field u Q (k,f) in the usual way, 

u Q (x, t) = Y, Mk, 0e ik x - (20) 

k 


Then the solenoidal Navier-Stokes equation for incompressible fluid motion takes the 
form (e.g. see McComb 1990) 



+ v 0 k 2 )u a (k,t) = \M a01 (k) J d 3 ju /3 (j,t)n 7 (k - j,t), 

(21) 

where 


M a0y (k) = {2i)- 1 {k 0 D ai {k) + fc 7 £> a/? (k)}, 

(22) 

and 



(23) 


A is a book-keeping parameter, which can be used to facilitate the iterative solution 
of the Navier-Stokes equation by expanding in powers of A. One then sets A = 1 at 
the end of the calculation. This procedure is well known in statistical and many-body 
physics and is usually referred to as the A-expansion. It is believed, however, that this 
expansion is divergent for the Navier-Stokes equation (Orszag 1977) nevertheless it 
is utilized in all RG theories. 

Turbulence almost invariably occurs because there is a mean rate of shear in a 
fluid. If we wish to study stationary, isotropic turbulence, then the artificial nature of 
this problem requires us to add hypothetical stirring forces to the right-hand side of 
the Navier-Stokes equations, in order to maintain the turbulence against viscous dis- 
sipation. Denoting the random force by / a (k, f), it is usual to specify its distribution 
as multivariate normal and to choose the force-force covariance to be 

</o( k, k, t ')} = D a/3 (k)W(k)S(t - t'), (24) 

where W(k) has dimensions of velocity 2 /time, and remains to be specified. Station- 
arity requires that the rate at which the stirring forces do work is the same as the 
rate of viscous dissipation e, or 



4.7rk 2 W(k)dk = 




(25) 


where 

roc 

e — J 2 u 0 k 2 E(k)dk. (26) 

It has always been normal practice in turbulence theory, to choose W(k) to be 
peaked near the origin, so that its arbitrary nature is only of importance at low 


12 



wavenumbers, and a universal energy spectrum can develop at large wavenumbers. 
More recently, the application of RG to stirred fluid motion (Forster et al., 1977; 
Yakhot and Orszag. 1986) has introduced theories which depend strongly on the 
choice of W(k). Accordingly, questions of how W(k) is chosen and justified must 
then be considered. In particular, in this context, the question of whether or not 
W(k) should depend on some characteristic length scale (such as the ultra-violet 
cutoff A) will surface later as a major issue. 

In order to try to carry out the RG algorithm given above in Section 2.2, we 
consider the Fourier components to be defined on the interval 0 < k < k 0 , where ko 
is the largest wavenumber present and is of the order of the Kolmogorov dissipation 
wavenumber (e.g. sec McComb 1990). Then, for some k\, such that ki < ko, we filter 
the velocity field at k = k\. Here, k\ = (1 — r))ko with 0 < rj < 1. This may be 
expressed in terms of the unit step functions, 


e~(k) 

e*(k) 


1 if 0 < k < k\ 

0 if k\ < k < k 0 ] 

0 if 0 < k < h 

1 if k\ < k < ko, 


allowing us to define the following useful filtered forms: 


(27) 

(28) 


u a (k,t)=$ ( k)u Q (k,t ), 
w t( k ’ 0 - 6 + {k)u a (k, t ), 

M-^(k) = r(*)M afl ,(k), 

Kf,( k) = 0*(k)M alh ( k). 

Then substitution of these forms into equation (24) allows us to decompose the Navier- 
Stokes equation into separate low- A; and high-/e forms, viz., 

(| + ^>-(k,t) = XM-^k) J d 3 j {up(i,t)u~(k- j,t) 

+2up(j,t)u+(k-j,t) 

+ttj(j,i)«+(k-j,i)}, 0 <k<h (29) 


r\ 

(^ + ^ 2 X( k ,*) = J d 3 j {up{i,t)u~(k- j,t) 

+2up(j,t)u+(\t- j,t) 

-Ht^(j,f)tt+(k - j,t)}. ki < k < k 0 (30) 

If we now try to carry out the first of the three steps outlined at the end of 
Section 2.2, for the lattice-spin Hamiltonian, we can adapt that procedure to the 
present (Navier-Stokes) case, as follows: 
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1. Solve equation (30) for u + . 

2. Substitute this solution into equation (29) for the u~ and do a partial average 
over the u + . 

3. Terms resulting from this procedure which are linear in u~ can be interpreted 
as contributing an increment to the turbulent viscosity. 

The main problems which arise due to mode- mode coupling, are easily understood: 

• First, the solution of (30) for u + contains terms in u~ . When these are sub- 
stituted into the right-hand side of (29), the result is a triple nonlinearity in 
the u~, which can then generate even higher-order nonlinearities in subsequent 
iterations. 

• Second, averaging out the high- wavenumber modes requires the hypothesis 

(■ u~u + u + ) = u~ (u + u + ), 

which cannot be strictly true, as u~ and u + arc just parts of the same velocity 
field and arc not statistically independent. < ... > represents the partial average 
over the high-wavenumber modes. 

It cannot be emphasized too strongly that these two problems arc fundamental 
stumbling blocks in applying RG to the Navier-Stokes equation. It is the way in which 
they are tackled which distinguishes one RG theory of turbulence from another. We 
will discuss in detail two methods for dealing with these problems: the recursive 
mode elimination and the conditional averaging method. In the approach to RG of 
Forster et al (1977), Fournier and Frisch (1993), which culminated in the theory of 
Yakhot and Orszag (1986), suitable assumptions axe made so that these problems are 
avoided. This approach is described in the Sec. 6 of this review. 

3.3 Iterative conditional averaging: a summary 

We now introduce a method of eliminating turbulent modes which is based on the 
use of a conditional average to distinguish between amplitude and phase correlation 
effects. It has its roots in the method of iterative averaging, which was developed 
over a number of years as a possible method of applying the renormalization group 
approach to real fluid turbulence (McComb 1982, 1986, 1990). However, an essential 
feature of the more recent work is the formal treatment of the conditional average and 
the development of methods of approximating its relationship to the usual ensemble 
average (McComb and Watt 1990, 1992; McComb, Roberts and Watt 1992). 

The basic idea is that the turbulent velocity field in wavenumber space may be 
decomposed into two distinct fields. One is a purely chaotic field ; while the other is a 
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correction field , and carries all the phase information. Application of this decomposi- 
tion to a thin shell of wavenumbers in the dissipation range allows the elimination of 
modes in that shell; with the usual mode-coupling problems being circumvented by the 
use of a conditional average. The (conditional) mean effect of the eliminated modes 
appears as an increment to the viscosity, with terms of order 77 2 being neglected, 
where rj is a dimensionless measure of bandwidth thickness, such that 0 < 77 < 1. 
An iteration (with appropriate rescaling) to successively lower shells, reaches a fixed 
point, corresponding to a renormalized turbulent viscosity. 

3.4 Recursive renormalization group: a summary 

If we return to the equilibrium Ising model, in Sec. 2.3, in which the Hamiltonian H 0 
involves only nearest-neighbor interactions, it was shown by Wilson (1975) that the 
effective Hamiltonians H n (n > 1) now involve both nearest-neighbor and diagonal 
nearest-neighbor interactions. Form invariance of H n is then imposed by assuming 
that the interactions of still higher order arc small. Similarly, in our form of recursive 
RG (Zhou ct ah, 1988, 1989; Zhou and Vahala, 1992, 1993a, b), as applied to the 
Navier-Stokes equations, form invariance can be imposed, when both the quadratic 
and the RG-induced cubic nonlincaritics arc retained, by truncating the A-cxpansion. 

The recursive RG theory extends Rose’s (1977) treatment of the linear problem of 
passive scalar convection to Navier-Stokes turbulence and eliminates the small scales 
recursively. Navier-Stokes equation include (1) a renormalized eddy viscosity and (2) 
a triple product in the fluid velocity. The eddy viscosity is calculated by means of a 
difference equation (Zhou et al. , 1988, 1989) or, in the limit 77 — » 0, by means of 
a differential equation (Zhou and Vahala, 1992). The resulting eddy viscosity term 
exhibits a mild cusp behavior in the renormalized momentum equation (Zhou et al., 
1988, 1989). The triple velocity products are absent only when there is a spectral 
gap between the subgrid and resolvable scales (Zhou ct al., 1988). Beyond the second 
iteration, the RG-induced triple products also contribute to the eddy viscosity in the 
renormalized Navier-Stokes equation. The recursivec RG, therefore, has potential for 
capturing a variety of dynamical features that depends on the interplay between local 
and nonlocal triad interactions (Zhou, 1993a, b; Zhou et al. , 1995). 

It is of particular interest to determine the effect of the new triple velocity product 
on the resolvable scale energy transfer and the difference, if any, from that of the 
usual Navier-Stokes quadratic velocity product, and we do this by considering their 
individual contributions to the eddy viscosity. Within recursive RG theory (Zhou and 
Vahala, 1993a), it has been shown that the triple velocity product in the renormalized 
momentum equation, w r hich produces a fourth order velocity product in the energy 
equation, removes energy from the resolved scales. The averaged fourth-order product 
can be decomposed into a product of averaged second-order products and formally 
takes the form of a gradient diffusion process with eddy viscosity v T . The triple 
term contribution to the eddy viscosity is zero as k — » 0 and increases rapidly as 
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k — > fc c . Here k c is the cutoff wavenumber which separates the resolvable and subgrid 
scales. The spectral eddy viscosity is simply the sum of the contributions from the 
momentum equation and the triple nonlinear term and it appears to be in qualitative 
agreement with the closure theory (Kraichnan, 1976; Chollet and Lesieur, 1981) and 
direct numerical measurement (Domaradzki et al., 1987; Lesieur and Rogallo, 1989; 
Zhou and Vahala, 1993a). In particular, it predicts the correct asymptotic behaviors 
of the eddy viscosity as k — ► 0 and k — * k c . The methodology has been applied to 
the passive scalar being advected by incompressible turbulence (Zhou and Vahala, 
1993b). 
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Part III 


Technical Aspects of selected RG 
approaches 

4 Iterative conditional averaging: details 

In this section we amplify our earlier summary of the method and show in some more 
detail how the iterative conditional averaging technique is used the implement the 
general strategy of the renormalization group for Navier-Stokes fluid turbulence. 

4.1 The conditional average 

First of all, we formulate the operation of taking a conditional average. Then we use 
this average in such a way that we can eliminate the triple nonlinearity referred to 
earlier. 

The idea is quite simple. Putting it at its most basic level, we select from the full 
ensemble of turbulence realizations a sub-ensemble, the members of which have their 
low-A; modes equal to u“(k, t ). Then we perform averages of functions of the u+(k, t ) 
over this sub-ensemble. 

However, there is more to it than this. If u(k, t) is the solution of the Navier- 
Stokes equation, corresponding to prescribed boundary conditions, then we arc faced 
with (in principle) a deterministic process and to prescribe u“ is to prescribe u + . 
That is to say, if u~ is invariant under our conditional average, then so also is u + . In 
order to get round this problem, we invoke the defining characteristic of deterministic 
chaos. This is to the effect that any uncertainty in the specification of the system 
will be amplified exponentially; so that as time goes on the difference between almost 
identical solutions will increase to the point of unpredictability. In the present case, 
we replace the concept of time going on, by the number of steps of the cascade 
in wavenumber. That is, our ideas about the turbulent cascade, and particularly 
ideas about localness of energy transfer, suggest that, if we prescribe conditions at 
wavenumber /c l7 then u^(k 0 ,t) will be unaffected, provided that A : 0 is much larger 
than k\. In other words, as the bandwidth becomes large (77 — > 1 ), the conditional 
average of becomes free of constraint and we can expect that 

{w + (k 0 ,i)) c -» ( u + (k Q ,t )). 

On the other hand, for ka — > k\, it is intuitively clear that the conditional average 
must tend to become effectively deterministic, with 

(u + (ko, t)) c -* u + {k 0 ,t). 
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Wc shall return to this at a later stage, but for the present it should be appreciated 
that the point at issue here is one of phase correlation. If two realizations are slightly 
out of phase at wavenumber k \ , then we can (given the nature of turbulence) expect 
that this phase difference will amplify exponentially, so that phase correlation will 
decline throughout the band as one moves from k\ to k 0 . Hence, for u + (k 0 , t) to be 
chaotic, despite the prescription u + (k!,f) = it - (ki,£), it is a requirement that the 
bandwidth be large enough. This imposes a lower bound on possible values of the 
bandwidth parameter r]. 

We may take account of these various aspects by introducing a fuzzy criterion 
for our conditional average. We now choose as our sub-ensemble, the subset of real- 
izations for which the low-A; modes differ from u~(k, t) by a small amount $~(k, t). 
Obviously d>~ is an arbitrary criterion (apart from the constraint that u~ + d>~ must 
be a possible solution of the Navier- Stokes equation) and should be chosen to satisfy 
the conditions which we wish to impose upon our average. 

Wc shall not go into this procedure any further here, and in order to simplify the 
algebra, we shall omit the from the equations. The interested reader who wishes 
to pursue the matter will find a full account in the papers by McComb and Watt 
(1992) and McComb et al (1992). 

4.2 Conditionally-averaged equations for high and low wavenum- 
bers 

Let us now denote the operation of taking a conditional average over the modes in 
the band k\ < k < ko by angle brackets, with a subscript 0. This notation per- 
mits the subsequent generalization to subscripts 1,2, ...,n, as we remove shells of 
wavenumbers progressively. Then, wc list the ideal defining properties of the condi- 
tional average as 

(u“(k, t)) 0 = u _ (k, <), (31) 

(u"(k, f)u“(k', t)) 0 = u“(k, i)u“(k', t), (32) 

and so on, for products of the low-wavenumber modes of any order. 

Wc now conditionally average both equations (29) and (30) which are, respectively, 
the low-A; and high-A; filtered Navier-Stokes equations. 

First, wc obtain the conditional-averaged NSE on the interval 0 < A; < k\, by 
averaging according to equations (31) and (32), to obtain 

+ i' 0 k 2 )u-(\t,t) = J d 3 j {up{j,t)u-(k-i,t) 

+2up{j,t)(u+(k - j,t)) 0 

+ (u^(j,A)a+(k- j,A))o}. (33) 
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Wc now repeat the steps just taken, but this time we apply them to equation (33) 
for the high-/;; modes. Thus, we get 

+ vok 2 )u+(k, t) = J d 3 j {up(],t)u-{k- j,f) 

+ 2 u^(j,/)u+(k - j,t) + Wg(j, t)u*(k — j, t)}, ( 34 ) 

and taking the conditional average of each term according to equations (32) and (33), 
gives 


(| + ^ 2 )(«^k,t))° = M+fa(k) J d 3 j {lip (j, /)u 7 (k — j, Z) 

+2u^(j,i)(u+(k-j,i))o 

+(up{j,t)u+(k -j,t)) 0 }. 

(35) 

4.3 First-shell elimination using the two-field decomposition 

Our objective now is to solve equation (35) for the conditional average { u p (j, t)u+( k — j, t)} 0 
and substitute the result back into (33), in order to have a closed equation for the low- 
wavenumber modes. In order to do this, we shall ultimately have to reckon with the 
need to relate conditional averages to full ensemble averages. Accordingly, we begin 
this section with the two-field decomposition which is our basis for this procedure. 

Let us write the exact decomposition: 

«J(M) = t£(M) + A+(k,t), (36) 

where u+(k, t) is any other realisation of our turbulent ensemble. In other words, 
u+(k, t) has exactly the same statistical properties as u+(k, t), but has no phase 
relationship to u“(k, t). It follows, by definition, that A+(k, t) is simply a measure 
of the phase difference (in the band of modes to be eliminated) between the two 
realizations. It also follows that, from the point of view of the realization under 
study, v+(k,t) is the purely chaotic part of the field and A+(k, t) is the correction 
field which carries all the phase information. 

Wc are now in a position to write down an expression relating the conditional 
average of the high-wavenumber part of the velocity field to its full ensemble aver- 
age. Taking the conditional average of both sides of equation (33), it may be shown 
McComb et al (1992), or it is intuitively obvious, that 

« ( k > 0)o = (Va (k, t)) + (A+ (k, t)) 0 . (37) 

Now wc seek a relationship between v + and u + , which is such that the conditionally 
averaged correction term (A + )o may be neglected as small. In other words, wc need 
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an ansatz for the correction term and naturally this will depend on the physical nature 
of the system that we are studying. 

In the case of macroscopic fluid turbulence, we axe guided by the well established 
idea that turbulent energy transfer in wavenumber takes the form of a cascade and is 
therefore to some extent local in wavenumber. In terms of our present approach, we 
take this to mean that, in any particular realisation, the effect on phase correlation 
of the mode-mode coupling is short-range in Ar-space. Thus, on a statistical picture, 
based on many such realisations, modes which are widely separated may be taken to 
be independent of each other. Hence, providing that the bandwidth parameter 77 is 
not too small, we can assume that u + (k 0 ,t) is independent of u + (ki,f), in the sense 
that we can write 

(u + (k 0 , t)) 0 = (u + (k 0 , t)) = <v + (k 0 , t)), (38) 

where the last step follows from the definition of v + , as another realization of the 
turbulence ensemble with the same statistical properties as u + , but with no phase 
relationship to u~ . 

This now leads us towards a natural ansatz for the relationship between v + and 
u + . Relying on the fact that we are dealing with a problem in continuum mechanics, 
we take v+ (k, t) to be given by a first order truncation of the expansion of u+(k,t) 
in Taylor series about k = k 0 , thus: 


t£(M) = Ma(k 0 ,t) + (k-ko).V fc u+(k,t) \ k=ko +0(r) 2 ). (39) 

Note that we conclude that terms of order rf have been neglected because the maxi- 
mum value of |k — k 0 | is 77 & 0 . Hence it follows that we have 

(A+(k,t)) 0 = 0(rj 2 ). (40) 


It also follows, therefore, that we are simultaneously imposing both upper and 
lower bounds on acceptable values of 77 . On the one hand, 77 must be large enough for 
us to assume that u + (ko,t) is independent of M + (k 1 ,t); while, on the other hand, 77 
must be small enough for us to neglect terms which are of order rj 2 in equation (39). 

Then, with all these points in mind, above equations yield for the viscosity acting 
on the explicit scales: 

V\ = v 0 + 8v o, (41) 

where the formula for the increment to viscosity is 


Sv 0 



3 L(k,j)Q+(|k-jl) 
J v 0 j 2 + i/ 0 |k - j | 2 


+ Ofa"*), 


(42) 


with 0 < k < ki,k\ < j, |k — j| < k () and Q+ is merely an extension of the spectral 
density to the v + field. The coefficient L(k, j) is given by 


L(k,j) = -2M p ^(k)M 0pS (j)D 6l (\k-j\) 

[t*(k 2 + f) -kj{ 1 + 2/x 2 )](l - n 2 )kj 
k 2 -I- j 2 — 2 kjn ’ 


(43) 
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where p is the cosine of the angle between the vectors k and j. 

We extend the procedure to further wavenumber shells, as follows: 

(a) Set u“(k, t ) = u Q (k, t) in equation (19), so that we now have a new NSE with 
effective viscosity i-'i ( k ) for Fourier modes on the interval 0 < k < k\. 

(b) Make the decomposition of (61), but this time at k = k 2 , such that w+(k, t ) is 
now defined in the band k 2 < k < k\. 

(c) Repeat the procedures used to eliminate the first shell of modes in order now 
to eliminate modes in the band k 2 < k < k\. 

In this way, we can progressively eliminate the effect of high wavenumbers in a 
series of bands k n + 1 < k < k n , where 


kn = (1 - v) n ko ; 0 < T) < 1, (44) 

with, by induction, the recursion relation for the effective viscosity given by 

iy n +i(k) = u n { k) + St/ n (k), (45) 


where the increment of order n takes the form 

, r , 3 . njoHgtojbfa + (i - + o(v 2 )} 

A> eJ 1 v„(j)P+v„( |k-j|)|k-j | 2 

Also, we may form an energy equation for the explicit scales, hence obtaining the 
renormalized dissipation relation, viz: 

[ kn 2 u n {k)E{k)dk = e, (47) 

Jo 

which may be compared with the unrenormalized where k n and v n should be replaced 
by k[) and u 0 , respectively. 

If we now assume that the energy spectrum in the band is given by a power law 
and make the scaling transformation 



k n +i = hk n , (48) 

where, for compactness, we define h = (1 — 77 ), it follows from equations (45) and (46) 
that the effective viscosity may be written 

u{k n k') = a 1 / 2 e 1 / 3 fc“ 4 / 3 i> Tl (A; / ) (49) 


where a is the constant of proportionality in the assumed spectrum. Now the recur- 
sion relation becomes 


with 


Z'n+l 


(k 1 ) = h i/z v n {hk') + h- A/z 5u n {k') 






'n(hj')j 12 + v n (hl')l' 2 


(50) 

(51) 
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( 52 ) 


for the wavenumber bands 0 < k! < 1; 1 < j\ l' < h~ l where /' = |k' — j'|, and 

Q' = h n/3 - _ A-i) + h.o.t. 

o 

Iteration of equations (50) and (51) reaches a fixed point and once this fixed point 
is found we can calculate the Kolmogorov constant by solving equations (47) and 
(49) simultaneously. One merit of taking a as a test is that it does have known 
experimental values, albeit scattered in the range 1.2 < a < 2.2. Our calculated value 
of the Kolmogorov spectral constant is a = 1.60 ± 0.01 independent of bandwidth in 
the range 0.25 < r/ < 0.45, and in good agreement with experiment. For values of rj 
outside this range, the calculated a diverges from the experimental value. At large 
values of 77, this is due to the breakdown of the first-order Taylor series approximation, 
while at small values, one is seeing the effects of mode coupling, which would invalidate 
the assumption that u(k 0 ) is independent of u(kx). 


5 Recursive RG 

5.1 Introduction 

We shall apply recursive renormalization group (RG) procedures to the problem of 
subgrid modeling. Subgrid modeling is necessary for the high-Reynolds number tur- 
bulent flows of interest because of the limitations of current and forseeable super- 
computers. A motivation is that the spectral transport coefficients (such as the eddy 
viscosity) determined from recursive RG theory can be compared to those arising 
from closure-based theories (Kraichnan, 1976; Chollet and Lesieur, 1981). It should 
be noted that the transport coefficients in these closure theories are determined over 
the entire resolvable scales, and are a function of k in the resolvable scales. 

In particular, we point out here that in e-RG, a small parameter e is introduced 
through the forcing correlation function. Yakhot & Orszag (1986) have to extrapo- 
late from e 1 to e — > 4 in order to reproduce the Kolmogorov energy spectrum. 
Furthermore, it is also necessary to take the distant interaction limit, k — ► 0. Thus, it 
is difficult to compare the wave- number dependent transport coefficients (Kraichnan, 
1976; Leslie and Quarini, 1979; Chollet and Lesieur, 1981), with that determined 
from e-RG. 

In recursive RG, no attempt is made to introduce a special form of overlapping 
as in conditional averaging, but one proceeds directly with standard averaging and 
handles the triple nonlinearity directly. The basic differences between the recursive 
and c RG procedures are that in recursive RG: 

(i) The e-expansion is not applied. 

(ii) The turbulent transport coefficients are determined for the whole resolvable 
wavenumber scales, 
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(iii) Higher order nonlinearities arc generated in the renormalized momentum 
equation and play a critical role in determining the transport coefficients. 

(iv) RG rescaling, as in conditional averaged RG, is performed. 

It should be emphasized that there arc two singular limits: h — ► 1 and k — » 0. A 
careful analysis must be done regarding these two limits and the associated averaging 
operations. We will address this issue here in the present review. 

5.2 Renormalized momentum equation 

For notation consistency, we shall use “+, — ” instead of “>, <” ■which arc typically 
employed in recursive RG. 

The first and third terms on the RHS of (29) arc symmetric in j and |k — j | in 
terms of their respective wavenumber constraints in wavenumbers. As a result, the 
distant interaction limit k — ► 0 has no effect on the existence of these terms, and 
these terms will give rise to the standard quadratic nonlinearity (first term of Eq. 29) 
and eddy viscosity (third term of Eq. 29). However, the second term on the RHS of 
(29) has the following constraint: j is in the subgrid while |k — j | is in the resolvable 
scales. Specifically, the consistency condition requires that, for small k, j satisfies 

j > k\ and j < k\ + kz 

where k • j = kjz. Since \z\ < 1, the range of integration must be O(k). 

Thus, the second term on the RHS of Eq. (29) can not contribute in the limit 
A; — > 0 since the integrand is bounded. Now it is well known that the higher order 
nonlinearitics are induced by this second term under discussion. Since this term is 
absent in the k — > 0 limit, we conclude that the higher order nonlinearies will not 
contribute to the renormalized momentum equations and recursion relation for the 
transport coefficients in the distant interaction limit , k — » 0. However, they will 
contribute to the renormalized Navicr-Stokes for 0 < k < hi. Atcr obtaining the final 
renormalized Navicr-Stokcs equation 


[d/dt + u(k)k 2 }u a (k, t) = f a ( k, t ) + M a01 { k) J d 3 ju 0 (j, t)u 7 ( k - j, t) 

+2 k) [ d 3 jd 3 j'-^M l3l3 , ll {j)u 0 >(y, t)u Y ( j - j', t)w 7 (k - j, t). (53) 

i/(k c )k c 3 J J 3 

by removing other subgrid shells iteratively, one can write down the corresponding 
wavenumber restrictions and perform the same analysis on the k — » 0 limit. 

5.3 Galilean invariance of the renormalized Navier-Stokes 

Here, we turn our attention to the question of the Galilean invariance of the renormal- 
ized Navier-Stokes equations. The importance of Galilean invariance in turbulence 
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modelling has been emphasized by Speziale (1985). To be consistent with the basic 
physics, it is required that the description of the turbulence be the same in all inertial 
frames of reference. The appearance of the triple nonlinear term, which is a function 
of the resolvable scales velocity fields, makes the property of the Galilean invariance 
of our recursive RG procedure not transparent. We now show that the renormalized 
Navier -Stokes equation is Galilean invariant (Zhou and Vahala, 1993b). 

The Galilean transformation is 


x — * x* — U 0 t* t —> t*. 
Here U 0 is an uniform velocity field. Thus, one has 


* T T v 

u = u - U 0 , tt = 


d_ 

dx 


d d_ _ d_ _d_ 
dx* dt dt* 0 dx* 


While the Galilean transformation for the Navier-Stokes equation in physical space 
is trivial, the Galilean transformation in wavenumber space is less obvious, due to 
the lack of differential operations. For convenience, we first review how Galilean 
invariance is preserved for the Navier-Stokes equation in the wavenumber space. 

Under the Galilean transformation, the LHS of the Navier-Stokes equation [cf. 
Eq. (21)] becomes 


dv* fk* t) 

1 + Uofiikfoi k*, t) + v 0 k* 2 [-U 0a 8(k*) + <(k*,f)] 
du* Ik* t) 

= + UvikfrX k*, t) + v 0 k* 2 <(k*,t) 

where in the last step, we have used the the 5 function property k* 2 d(k.*) = 0. 

Also, under the Galilean transformation, the RHS of the Navier-Stokes equation 
becomes 


M afil {k*) J d 3 j[u* 0 u*,t) - iM(j*)K(k* - j*.0 - ^^k* - j*)] 

= M a 0 y (k*) J d 3 j*u* 0 (}* ,t)u*(k* — j*,f) + iUopkpUaCk* ,t) 

where we have used the property of the 6 function, the incompressible condition, and 

M a01 (k*)u O0U ;( k*, t) - u O0 k; u * a (k*, t)/(2i) 

Thus, as expected, the Navier-Stokes equation is invariant under a Galilean transfor- 
mation due to the cancellation of the second term on the RHS. 

To show that the renormalized Navier-Stokes equation is invariant under a Galilean 
transformation, we need only consider the recursive RG induced triple nonlinear term, 
denoted by N St- 
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NS t = 2 M a 0 y (k) J d 3 jd 3 j'(i-) 4 / 3 ~ u 0 ' (j - j', t)uy( j', t)tt 7 (k - j, t) (54) 

It is important to note that j is in the subgrid. 

Under a Galilean transformation, Eq. (54) becomes 


JVSJ. = 2 J d 3 i'rfV(f ) 4/3 MJ* - jV) " V v S(s‘-i 


(K)j' 


K-av) - iv<ariK(k - - r,t) - tw - m 

Since j* is in the subgrid scale, while f* and k * are in the supergrid, 5(k* — j*) 
and 5(j* — j'*) can never be simultaneously satisfied. As a result, 






(55) 

Now only one term in Eq. (55) could violate the Galilean invariance of the renormal- 
ized Navier-Stokes equation. However, 


-r.i) 


j d 3 fs(Du;,( j* -r,o u- s ,(y,t) 

This is not permissablc since u 0 , = u—* 0 > and j* is restricted to the subgrid. Thus 
NSt — NS^- Hence the triple term is Galilean invariant. 


5.4 The effect of cubic nonlinearity 

We consider the contribution of the triple nonlinear term in the renormalized mo- 
mentum equation to the eddy viscosity. The second moment for the velocity field is 
defined as 


U a0 { k, t) =< u a ( k, t)u 0 (- k, t ) > . (56) 

The time evolution of U a p(k,t) is 

aC/ °^ k,f) = -2u(k)k 2 U^(k, t) + 1%, (k, t) + T^( k, ()• (57) 

In this equation, T® a (k,t) is the standard energy transfer from the quadratic nonlin- 
earity. In contrast, Tj Q ( k, t) = — 2u T {k)k 2 E{k) is the energy transfer arising from the 
RG induced triple nonlinearity. It is readily shown that 
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v T {k) - 


(58) 


7£(M) = _J_1 /*+*»... 9)1^ — jl~ y - 2 j v+1/3 

2E(k)k 2 2u(k c ) k 2 Jk c 3 Z i/(k — j) 

It has been shown that that Vr(k) is the major contributor to the cusp-like be- 
havior of the spectral eddy viscosity as k — ► 0 (Zhou and Vahala, 1993a). This type 
of term was found to be the major contribution to the strong cusp in the spectral 
eddy viscosity found from the closure models (Kraichnan, 1976; Leslie and Quarini, 
1979; Chollet and Lesieur, 1981). 

Rose (1977) discussed the role of the triple nonlinear terms in physical space. He 
pointed out that it represents the possibility of an exchange of scalar eddies between 
the resolvable and subgrid scales. This effect is an inherent property of measurements 
made on the passive scalar system with instruments which have a spatial resolution 
limited to an eddy width size greater than 1 jk c . 


5.5 Difference equations for the renormalized eddy viscosity 

After the removal of the (n + l) lh subgrid shell, the spectral eddy viscosity in the 
renormalized momentum equation is determined by the recursion relation 


where 


and 


^n+l(^) ^ / n(^*') T ^^n(^) 


5v n {k) = 


Do^ f ,3- LjkEhQ) jk~j| y _ 

k 2 hJ J "h(j)j 2 v( |k-j|)|k-j| 2 ’ 


(59) 

(60) 


L(k,j,q ) = ~ 


kjjl - z 2 )[ z < 7 2 

Q 2 



(61) 


with k • j = kjz and q = |k — j|.. This difference equation, after rescaling, has been 
solved by Zhou ct al. (1988, 1989) and fixed points were readily determined for 
finite h < 0.7. Note that the spectral eddy viscosity shows a mild cusp a s k —> k c , 
in qualitative agreement that that of closure theory. However, it was very difficult 
to determine fixed points for finer subgrid partition factor h > 0.7. In the next 
subsection we shall pass to the differential subgrid limit h —> 1 and determin an 
ordinary differential equation (o.d.e) for the renormalized eddy viscosity over the 
entire resolvable scale which can be readily integrated. 


5.6 Differential equations for the renormalized eddy viscosity 

The differential limit, h — ► 1, is a singular one and this point has been discussed 
recently (Zhou and Vahala, 1992). In particular, it is related to the assumption of 
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local versus non-local interactions in k. In this section we will calculate the eddy 
viscosity under the differential equation limit for recursive RG. 

For recursive RG we will find that the differential equations hold throughout the 
resolvable wavenumber range 0 < k < k c . This should be contrasted with the e - RG 
eddy viscosity differential equation which is valid only in the k — » 0 limit. 

5.6.1 The differential equation limit, h — > 1 

We now derive the differential equation from which the transport coefficients for finite 
k, 0 < k < 1, arc determined. The o.d.e. in the distant interaction limit (k = 0) 
will be derived in the next subsection. After the rescaling, we rewrite the recursion 
relation in the form 


u n+ i(k) - h iy+1)/3 u n (hk) = h (y+l)/3 5v n {hk). (62) 

For h — * 1, the number of interation n — > oo. Similarity consideration leads to 

v n +i(k) — * is(k), n — > oo. (63) 

The LHS of Eq. (62) becomes 

u{k) - [1 - r)} iy+1)/3 v{k{ 1 - 77 )] -» v[-^- + yyKO + Ofa)] (64) 

As noted earlier, the partial average of Rose must be employed in order to ensure 
the existence of the differential limit. The partial average is introduced since the 
distinction between the resolvable and subgrid scales becomes fuzzy in the limit of a 
differential subgrid partitioning, h — ► 1 . 

Following Rose, we first change the variable from j,z to j,q , with djdz = 
( q/kj)djdq , so that the RHS of Eq. (62) becomes 


5v{k) = J djdq(jt) 
J djdq(-jp) 


L(k, j, g) 

kj'v{j)k 2 v (|k — j|)|k— j| 2 |k — j|^ 


L(k,j, q) 


kj v(j)k 2 iy(\k — j|)|k — j| 2 |k — j|y k, 




V 


l 


L(k,l,q) f 
dq , ^ — 9 / 1 \ _ 1 + V I 

J 1 <7 


i< 9 <i-i-fc k 3 u 2 (l)q y 


<j<\+k 


j.L(k,j, 1) („_2)/3 

T J 


where 


L(k, 1, q) = k(l - z 2 )(k - zq 2 )/q 2 

and 

L{k,j, 1) = kj( 1 - z 2 ){kj - z). 


(65) 


(66) 

(67) 
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As a result, the fixed point renormalized eddy viscosity u(k) is determined from 
the o.d.e. at 0(A) 

+ B„«] (68) 

where 



LjK hq) 

qy-i 


(69) 


1 /* 1 -f- k 

Bv{k) = Vh djL{k ' j ' 1)j ' (v " 2)/3 (™) 

Here, z is evaluated at j = 1 and q = 1, respectively in the L(k,j,q) expression. The 
o.d.e. for the momentum equation eddy viscosity is readily solved (Zhou and Vahala, 
1993a). 


5.6.2 Differential equations in the k — » 0 limit 

In the limit k — > 0, we have seen that the triple nonlinearities induced by RG do not 
contribute to the eddy viscosity. As a result, the recursion relation will now contain 
only the usual quadratic contribution. We further simplify the analysis by taking 
the standard subgrid linear propagator Gj^flk — j|) = [J^ + ^hdk — j|)] <V(|j|) as 

k — > 0. 

The limits of the integration are given by 


1 < j - kz < 1 //, 1 + kz < j 

< 1 / f + kz 

(71) 

Thus, one has 




8» n (k) = S- 

-F-G 


( 72 ) 

where the integral limits for these terms 

are 



r'/f ri 

h di L dz 

for 

S 

( 73 ) 

r 1 rl+kz 

I dz L dj 

for 

F 

( 74 ) 

f0 rl/f 

/ dz dj 

J - 1 Jl/f+kz 

for 

G. 

(75) 


Terms F and G are the corrections to the symmetric term S. They are important 
for a finite bandwidth /. However, it is easy to show that F + G = 0 for / — ► 1 in 
the k — * 0 limit. Hence 


s ^ k) ~ 


15 v 2 {l) 
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while the LHS of the same equation yields 


. If* 


y + 1 


v(k), 


as 


since is bounded as k 
Thus, as k — + 0, 


0 . 


u(k —* 0) 


3 8 1 

y + 1 15 i/ 2 (l) 


(77) 


Again, the o.d.e. for the momentum equation eddy viscosity is readily solved. We 
observe that the eddy viscosity has a similar plateau structure as k — ► 0. As k — ► k c , 
eddy viscosity displays a weak cusp like behavior as k — > k c . In this case, those curves 
arc similar to that of Zhou et al (1988; 1989). 


5.7 Recursive RG Spectral eddy viscosity 

The spectral eddy viscosity is simply the sum of the contributions from the momentum 
equation and that of the effect of the RG induced triple nonlinear term in the energy 
equation. It is apparent that our calculation is in qualitative agreement with that 
from the closure theory (Kraichnan, 1976; Chollct and Lesicur, 1981). and direct 
numerical measurements (Damaradski et al., 1987; Lesieur and Rogallo, 1991; Zhou 
and Vahala, 1993). In particular, it predicts the correct asymptotic behaviors of the 
eddy viscosity as k — » 0 and k —> k c (Kraichnan, 1976). 


5.8 Numerically evaluated eddy viscosity 

These conclusions can be tested directly using numerical simulation databases. In- 
deed, energy transfer and spectral eddy viscosity can be analyzed using results from 
direct numerical simulations by introducing an artificial cut at a wavenumber k c that 
is smaller than the maximum resolved wavenumber k m of the simulation. With this 
fictitious separation between the subgrid and resolvable scales, it is possible to eval- 
uate the effect of the subgrid k c < k < k rn on the resolved scales k < k c . We form an 
energy equation from the momentum equation and introduce the following notation: 
T >K (k ) and T >> (k.) represent the spectrum of energy transfer to mode k resulting 
from interactions with one and both modes above the cutoff k c respectively. Measure- 
ments of numerical simulation databases indicate the following (Zhou and Vahala, 
1993a): 

• T >> (k) removes energy throughout the resolvable scales in a manner consistent 
with the notion of eddy viscosity. 

• T ><: (k ) removes energy from the last resolved octave that was transferred there 
by the resolved scale transfer; that is, it allows the local flow of energy through k c . 
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It is the most important subgrid effect near k c and accounts for most of the energy 
flow from the resolved scales. 

The subgrid spectral eddy viscosity r' >> (fc) and u >< (k) can be determined from 
T >:> (k) and T ><: (/t) for a given energy spectrum, E(k). Specifically, u >:> (k) = 
—T >> (k)/2k 2 E(k) and u >< {k) = —T >< (k)/2k 2 E(k). Two important features of 
the quadratic contribution u >>r (k) should be stressed. First, its positive constant 
asymptote at small k indicates that the concept of modeling this contribution as an 
eddy viscosity in analogy to the molecular viscosity is plausible, and second, its value 
decreases monotonically as k increases toward k c . This indicates that if we include 
only the contribution of quadratic velocity products, there is no eddy viscosity cusp at 
the cutoff k c . The most important feature of i/ >< (fc) is the sharp increase at k — ► k c . 

6 Reconsideration of the YO theory 

Robert Rubinstein 

Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, VA 23681 


The first comprehensive application of RG methods to turbulence was the Yakhot- 
Orszag (1986) theory (YO). This theory built upon earlier work by Forster et al 
(1977) and Fournier and Frisch (1983), but completed the development suggested by 
its predecessors with impressive successes including simple analytical calculations of 
the Kolmogorov constant 1 and several constants of interest in turbulence modeling 
(Speziale, 1991). Nevertheless, this theory continues to generate controversy. Three 
points will be addressed in this section about which particularly strong objections 
have been raised: 

1. YO’s treatment of mode elimination 

2. the “correspondence principle” 

3. the e-expansion and distant interaction limit 

This section draws heavily on work of Woodruff (1992, 1993) which emphasizes the 
connections between the YO theory and Kraichnan’s (1959) direct interaction ap- 
proximation (DIA). Although this viewpoint tends to de-emphasize the importance 
of mode elimination characteristic of renormalization group theories, it is consistent 
with the remark of Eyink (1994) that YO is not a true RG theory in any case. 

In this Section, the following notation will be used: 

k = (k, f!) k = (p, u) q = (q, Q — u) 

3 Thc Kolmogorov constant had been computed earlier by Kraichnan (1964) by numerical inte- 
gration of the LHDIA closure. 
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where the triangle condition 


k = p + q 


holds throughout. 

6.1 YO’s treatment of mode elimination 

It was noted at the end of Sect. 3.2 that iterated mode elimination leads to two 
fundamental problems: it generates nonlinearities of higher order than appear in 
the Navier-Stokes equations, and it requires some simplifying assumptions to permit 
averaging over the high wavenumber components of the motion. These difficulties 
arise because, by itself, mode elimination, like the field theoretic functional formula- 
tions of the turbulence problem, does nothing more than reformulate the equations 
of motion. In order to make progress, statistical hypotheses must be introduced. 
This circumstance is not surprising: for example, equilibrium statistical mechanics 
requires Gibbs’ hypothesis, a very strong assumption which certainly does not follow 
from Newtonian mechanics. 

Nevertheless, subsequent work has shown that the original treatment of mode 
elimination by YO was not entirely satisfactory. For example, their proposal that 
higher order nonlincarities arc negligible in a certain perturbative sense, was contra- 
dicted by Eyink (1994). 4 But by placing the YO theory in the setting of DIA-likc 
closures, these problems are understood by invoking the statistical hypotheses of these 
closures. Thus, higher order nonlinearities can indeed be considered, but in the con- 
text of higher order versions of DIA (Martin et al , 1978). Accordingly, the absence 
of these nonlinearities in the YO theory merely reflects the order of approximation 
chosen, and requires no further justification. 

The treatment of averaging over the small scale motion, which appears to proceed 
as if motions of different scale were independent, has been discussed in Kraichnan’s 
(1959) original presentation of DIA: the same requirement arises in the derivation of 
the DIA response equation. Hcuristically, the DIA closure assumes that the velocity 
field is only weakly non-Gaussian. Therefore, the motions of different scales, which 
arc uncorrelated because of the kinematic hypothesis of homogeneity, are independent 
to leading order. This statistical hypothesis also justifies the breakup of fourth order 
moments into products of second order moments. A theory in which the motions 
of different scale are not independent to leading order, but in which fourth order 
moments arc treated by the quasi-Gaussian hypothesis, must be carefully formulated 
to avoid inconsistency. 


4 YO had asserted that higher order nonlinearities are ‘irrelevant’ as this is understood in Wilson’s 
(1974) theory when e — 0, although they are only marginal when e — 4. Eyink (1994) demonstrated 
that, on the contrary, higher nonlinearities arc marginal regardless of e. 


31 



6.2 The correspondence principle 


The starting point of YO’s analysis is the Navier-Stokes equations driven by a random 
force: 

-iSlUi(k) - M imn (k) [ dpdq u m (p)u n {q) = fi(k) (78) 

Jk=p+q 

where the Gaussian random force f t is characterized by its correlation function 

< fi(k)fj{k') >= 2D(2n) d+1 k- y D lj (k)5(k + k') (79) 


This force is white noise in time. The exponent y is treated as a variable for purposes 
of the subsequent e-expansion, in which e = 4 + y — d and d — 3 is the dimension of 
space. The analysis leads to a Kolmogorov spectrum when y = 3 or e = 4; this is 
therefore the case of physical interest. YO’s conclusion is that the nonlinear term the 

Navier-Stokes equation is replaced, in the limit of infinite Reynolds number, by the 
combination of random forcing by fi and a scale dependent viscosity v(k), so that 


—iQui(k) + u(k)k 2 Ui(k) = fi(k) 

Mode elimination is used to obtain the recurrence relation 

dv(k) D 

dk u 2 k 5 


(80) 


(81) 


where the constant A is computed from the theory. 

Whereas it is generally agreed that Eq. (78) provides a plausible model of isotropic 
turbulence provided the random force / is concentrated at large scales and therefore 
provides an energy source, the introduction of a force acting on all inertial range 
scales appears to lack fundamental justification. YO’s model Eqs. (79) and (80) can 
be compared to the DIA Langevin model (Kraichnan, 1976) 


— iQu(k) + T](k)ui(k) = fi(k) 


(82) 


where the damping function and force correlation are expressed in terms of the DIA 
response and correlation functions G and Q by 

r](k) = 2iM rmn (k) [ dpdq D mr {p)D ns (<i)G(p)Q{q) (83) 

Jh=p+q 

F (k) = < fiityfrik 1 ) > /<5(k+k') - -4M imn (k)M jrs (k) [ dpdq D mr (p)D ns (q)Q(p)Q(q) 

Jk=p+q 

(84) 

Eq. (82) is a generic model in statistical mechanics which replaces the effects of an 
infinity of nonlinear interactions on any one mode by a random force acting against 
a generalized damping; DIA applies this description to a problem which is far from 
thermal equlibrium. The “fixed point” RG model Eqs. (79) and (80) formally resem- 
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bles the DIA Langevin equation model. But the damping i/(k)k 2 in in the RG model 
Eq. (79) is Markovian, so that 

rj{k) = r)(k) (85) 

only and the forcing in Eq. (80) is white noise in time, so that 

F(jc) = F(k) (86) 

only. Neither of these conditions holds for the DIA Langevin model. 

To investigate the connection between these models, write following Kraichnan 

Q(p) = Q{p)R(p) (87) 

where R is the time correlation function. Perform the frequency integration in Eq. 
(84) and evaluate the result in the long time limit in which Q = 0. This limit 
corresponds to observing the system over times long compared to any characteristic 
correlation time of the true DIA random force. The result is 

F(k) = -4M imn (k)M jrs (k) f dpdq D mr (p)D ns (q)Q(p)Q(q)Q(k,p, q) (88) 

^k=p-fq 

where 

/ oo 

duR(p)R(q) | n =o (89) 

-OO 

In this limit, the random force is white in time. Further, in Kolmogorov scaling, 

Q{\p) = X~ n/3 Q(p) (90) 

e{Xk,Xp,\q) = X- 2 ' 3 e{k,p,q) (91) 

consequently, the scaling dimension of the random force is found to be -3: 

F(Xk) = X ~ 3 k (92) 

Formally, in the long time limit, the random force in the DIA Langevin model has 
the same space-time correlation as the force postulated at the outset by YO. 

It should be noted that the power counting which leads to Eq. (92) is purely 
formal, since the actual force correlation integral in Eq. (84), like the integral of Eq. 
(83), is infrared divergent when evaluated for an infinite Kolmogorov inertial range. 
We recall that these divergences actually cancel in the DIA energy equation; however, 
the assumption of a -3 force in DIA requires a priori infrared regularization. 

That the -3 force is natural in the context of any steady state far from equilibrium 
with a constant flux of some invisicid invariant is also suggested by the derivation 
(Rubinstein, 1994a) of Bolgiano scaling inertial range for buoyant turbulence by ap- 
plying the YO formalism with forcing of the temperature equation only. From this 
point of view, the introduction by Lam (1992) of a distinguished infrared scale in the 
RG force is debatable: it corresponds to a loss of the locality of the inertial range 
postulated by Kolmogorov. 
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6.3 The e-expansion and the distant interaction approxima- 
tion 

The e-expansion is the subject of an especially large number of re-evaluations and re- 
considerations, among which are Ronis (1987), Lam (1992) and Wang and Wu (1993). 
In YO’s original presentation, the e-expansion is an expansion about a logarithmically 
divergent theory. An interesting alternative was suggested by Carati (1990a), who 
suggested expanding about a theory with vanishing energy transfer (Fournier and 
Frisch, 1978). Here, this expansion will be considered, following Woodruff (1992), as 
an approximation in DIA. 

To complete the transition from the DIA Langevin model to the YO theory, further 
approximations are required. They are 

• (a) evaluate the DIA integrals in the distant interaction limit in which k/p, k/q —> 

0 

• (b) Markovianize the damping 

• (c) introduce an infrared cutoff so that the integrals in Eqs. (83), (84) are 
restricted to p > k and q > k only 

• (d) evaluate the amplitudes using the e-expansion 


It has been emphasized by Woodruff that these approximations are closely related. 
First, as noted by Kraichnan (1987), the e-expansion is an expansion about a theory 
in which distant interactions are dominant; accepting this point provisionally, we 
outline how the distant interaction limit brings about the Markovianization of the 
damping and forcing. 

Let the real function i/(£),0<£<oo satisfy 

roc 

H{ 0) = 1 , H{0 < 1 for £ > 0, / H{Z)d£ < oo (93) 

Jo 

Then standard properties of delta functions imply 

XH(X(t — s )) ~ S(t — s) for A — > oo (94) 


Rewrite Eq. (83) in the time domain, and evaluate the wavevector integrals in the 
distant interaction approximation in which k — ■* 0, p, q — > oo. Then 

r?(k, t, s) = / dpdq B( k, p, q)G(p, t, s)Q(q, t, s ) 

Vk=p+q 


J dp {k m - — (k, p, p)G{p, t, s)Q(p, t, s) 
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-£(k, p, p )G(p, t, s)k m p m p 1 ^(p, t, s)} 

where Z3(k, p, q) denotes the product of projection operators in Eq. (83). Assuming 
time stationary similarity forms 

G{p, t, s) = G(p r (t-s )) 

Q(p, t, s ) = i?(p r (f - s))Q(p) 

the properties Eqs. (93) of H may reasonably be postulated of the product GR. 
Therefore, Eq. (94) implies that in this limit the damping is Markovian 

*/(k, t,s) = 5{t-s)rj( k) 

and the DIA response equation implies that the Green’s function is exponential, 

G(k, t, s ) = exp [(s — t)/ 7 (k)] for t > s (95) 

Likewise evaluating the force correlation Eq. (84) in the distant interaction limit 
implies that the forcing is white noise in time: 

< >= S(t - s)<5(k + k')Fjj(k) (96) 

Computing the two-time correlation function from the relation 

Qij( k, t, s)6(k + k') 

= [ dri G(k, t, r j) f dr 2 G{k',t, r 2 )x < f i (k,r 1 )f j (k\r 2 ) > 

Jo Jo 

using Eqs. (95), (96) shows that the fluctuation dissipation relation 

Q(k, t, s) = Q(k)[G(k, t, s) + G(k, s, t)] (97) 

expressing the time dependence of the correlation functions in terms of the response 
function is also valid in this limit. 

These simplifications of DIA permit analytical evaluation of the inertial range con- 
stants. Although values of these constants could be inferred from numerical solutions 
of DIA, say for decaying turbulence, it is natural to attempt analytical evaluation as 
well. 

Introducing Eqs. (95) and (96) with the Kolmogorov scaling forms 

E(p) = C K e 2/3 (fc“ s/3 (98) 

Mp) = Coe' V /3 

into the DIA response equation integrated over all time separations, 

C 2 f n -11 / 3 

cic = L +tl dPdq 2 * Mrmn ( k ) D rnr ( p ) D ns ( q ) ^ 2/3 — ^ 2/3 ^2/3 
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Integrating the single time equation for the correlation function with respect to 
wavenumber k leads to a second equation, which for Kolmogorov scaling gives 


Cd 

Cl 


= .1904 


( 100 ) 


This method of evaluating the inertial range constants Cd and C*- fails because of 
the well-known divergence of Eq. (99) at low wavenumbers. The simplest infrared 
regularization which is consistent with Kolmogorov scaling is to restrict the region 
of integration to triads satisfying p > ak. Values of Co(ct) and Ck{ol) have been 
tabulated by Leslie (1972). 

The e-expansion can be considered as a method of infrared regularization by an- 
alytic continuation. Namely, continue to assume Eq. (97), but replace Eq. (98) by 
the general form 

E{p) = CxD^k 1 - 2 *' 3 (101) 

The scale independence of the integrated response equation demands 


flip) = CoD^k 2 -^ 3 


( 102 ) 


The units of D, consistent with Eq. (79), make these equations dimensionally correct. 
Substituting these scalings in the integrated response equation gives the e-dependent 
form of Eq. (99), 


Cl 

Ck 


I dpdc\ 2iM rrnn (k)D rnr (p)D ns (q) 

/k=p-bq 


V 


-1 — 2e/3 


(p2-e/3 + gr2— e/3^2 — e/3 


(103) 


The integral in Eq. (103) is ultraviolet divergent when e < 0 and is logarithmic when 
e — 0. Woodruff observes that it is reasonable to evaluate Eq. (103) for e > 0 by 
asymptotic expansion about e = 0. This expansion greatly simplifies the integration 
since the ultraviolet divergence for e < 0 implies that the integral is dominated by 
distant interactions, namely by wavcvector triangles such that p.q —> oo. In this 
limit, a simple analytical evaluation of the integrals is possible. The calculation gives 

-pp- = — ^4(c) = f Aq + A\t + • • • (104) 

Uk e e 


where 



The constant /t_i is distinguished since it is the only one in the scries Eq. (104) which 
has been evaluated exactly in two senses. First, increasing the number of “loops,” 
that is, considering terms in the perturbative solution of the Navier-Stokes equations 
with a larger number of force correlations, will correct A p only for p > 0. It can also 
be shown (Rubinstein, 1994b) that even at the one loop level, correcting the distant 
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interaction approximation by power series expansions in k/p also only corrects A p for 
p > 0. Accordingly, it is reasonable to evaluate Eq. (104) by taking the leading term 
only. Setting e = 4, 

C| = _3_ 

C K 20 

which is easily shown to be equivalent to YO’s calculation. 

It is sometimes claimed that YO evaluate the inertial range constants by setting 
e to four and zero at different places in the same equation. Therefore, it must be 
emphasized that in this calculation, e is never set to any value but 4 - The analytical 
procedure which leads to Eq. (104) is entirely routine: it is the evaluation of the 
leading term in an asymptotic expansion, not a novel procedure unique to YO. 

The e-expansion was described earlier as an infrared regularization necessary to 
evaluate the right side of Eq. (99), which diverges when p — * 0. Triads with p ~ 0, q ~ 
k correspond to sweeping of modes with wavevector Ik!--/,' by modes of much larger 
scale. The dynamic significance of this divergence has been elucidated by Kraichnan 
(1982). This divergence is removed in YO, and indeed in all renormalization group 
approaches by focusing exclusively on interactions for which p : q > k. 

In fact, the e expansion is constructed so that when e = 0, the dominant inter- 
actions actually are the distant interactions for which p, q — + oo: in this case, the 
integral in Eq. (99) is logarithmically divergent in this limit. However, as Woodruff 
(1993) notes, the integral becomes infrared divergent when e = 3, and the analytic 
continuation from e = 0 to e = 4 in the YO theory becomes problematic. Thus, al- 
though it is satisfying to be able to compute the inertial range constants, and even to 
obtain satisfactory values by a straightforward computation, the fact remains that the 
analytic continuation which underlies the calculation requires justification. Moreover, 
Woodruff also suggests that one might attempt an e expansion about this infrared 
divergence. Not unexpectedly, the results are quantitatively unsatisfactory, but this 
possibility suggests that the expansion about e = 0 is not the only one possible. 

Another objection to this procedure can be raised in connection with Eq. (100): 
the constant has been obtained by exact evaluation of the triangle integrals making 
neither the e-expansion nor the distant interaction approximation. However, the 
integral can be shown to be ultraviolet divergent for e < 4 and logarithmic exactly 
when e = 4. Thus, there is no possibility of an e expansion for this integral, which 
must be evaluated exactly. 
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7 Appendix: Brief Description of RG Publications 
Relevant to YO 

7.1 The work of Kraichnan 

Kraichnan (1987) provides an extensive discussion on YO. First, he pointed out the 
logical distinction between two procedures: the multiple scale elimination and the 
e-expansion. First, the evaluation of eddy viscosity by perturbative elimination of 
successive small spherical shells of high- wavenumber modes can be carried out by 
specifying an actual energy spectrum (Rose, 1977) rather than introducing a forcing 
spectrum. Kraichnan discussed next the method of “e-expansion” . In this procedure, 
the properties of a E(k) ~ /c 1_2t/3 spectrum range are examined through an expansion, 
in powers of e, about the properties of a spectrum E(k) ~ k. 

Kraichnan’s analysis is based on estimation of the qualitative nature of eddy 
damping in a spectrum of the form 

E(k) oc k s (k c < k < kd ). (105) 

The corresponding eddy damping is given by RG as 

v{k\p) oc ( p/k )~ e / 3 (0 < e < 4), (106) 

where s = 1 — e. For the Kolmogorov inertial range spectrum, e = 4, the eddy 
damping is local for any positive e since v(k\p ) — >• 0 as p — ► oo. In YO, the explicit 
calculation of eddy damping is first made in the near neighborhood of this reference 
spectrum. Now distant interactions really are weakly dominant for E ( k ) ~ k. At 
s = 1 or e = 0, Eq. (106) is replaced by 

u(k\p) oc ln(A; c /p). (107) 

The results are then mapped to the Kolmogorov spectrum by taking e — ► 4. In 
the meantime, YO retained terms through first order in e, whereas 4 is not a small 
number. Hence, the convergence properties of the expansion are unclear. 

Kraichnan (1987) also stressed that this type of estimation can be done completely 
by dimensional analysis (Fourier and Frisch, 1978; 1983; Kraichnan, 1982). The only 
element of the Navicr-Stokes equation involved in these discussions arc the overall 
energy conservation by nonlinear terms and the coefficients of interaction of distant 
wave vector triads. 

7.2 The work of Teodorovich and Wang &; Wu 

Teodorovich (1987, 93, 94) used the field-theoretical method to re-evaluate the results 
of YO. We will discuss two issues discussed in his most recent publication in 1994. 
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The first issue investigated by Teodorovich is a reassessment of the eddy viscosity 
calculation of YO. This is exactly the same subject examined in Wang & Wu (1993). 
Recall that in YO the parameter e is used as an small pertuterbation expansion and 
is an important parameter in YO’s eddy viscosity expression. YO then took the 
limit 6 — ► 0 when they compute the coefficient but set e == 4 in the powerlaw expo- 
nent. However, both Teodorovich and Wang & Wu (among many others workers) 
noted that there is no mathematical justification for doing this, and the choice of 
e = 4 leads to unacceptable results. Furthermore, they found an algebraic error in 
YO. Specifically, in the integration over wavevector q, YO introduced a substitution 
q — » k + q/2. However, YO did not take this substitution into account in the cor- 
responding transformation in the associated domain of integration. Correcting this 
algebraic error, they found an eddy viscosity which is independent of the parameter 
e. Hence, the eddy viscosity in Teodorovich and Wang & Wu is independent of e, 
but has the numerical value of the eddy viscosity as YO who took the e — > 0 limit. 
Without e dependence in the eddy viscosity, Teodorovich and Wang & Wu propose 
that a good agreement with experimental data can be obtained. 

It is important to note the limit e — ■> 0 in YO was not only used to obtain an 
acceptable value for the eddy viscosity, but the limit e — > 0 also plays an important role 
in eliminating higher order nonlinearities in the RG procedure. However, Teodorovich 
and Wang & Wu do not address the closure problem if e = 4 is maintained throughout 
the analysis. 

The discussion now leads naturally to the question on how the local and nonlocal 
interactions are being treated by YO a second major issue in the work of Teodor- 
ovich. For a forcing correlation which is a function of e, an energy spectrum can be 
deduced with an e dependence. As discussed already, Kraichnan (1987) found that 
the interactions are local when e — > 0 but nonlocal when e = 4. Hence the results 
computed at one limit may not be used at another limit. Teodorovich claims that 
e — > 0 corresponds to the dominance of local interaction. An analytic continuation 
in e from a pole-type singularity at e = 0 to the point e = 4 means neglecting the 
nonlocal interactions with the large-scale modes (Teodorovich, 1994). His analysis is 
not transparent. 

7.3 The work by Ronis 

Ronis (1987) analyzed a model for randomly stirred fluids using RG on a path-integral 
representation of the Navier-Stokes equation. Unlike YO, he found that a choice of 
the random force correlation exponent (y = —1.5851 in three dimensions) is needed 
to give the Kolmogorov 5/3 law at high wavenumber. 

Of particular interest, Yuan and Ronis (1992) discussed the issue on the generation 
of random force. They noted that YO and Ronis (1987) ignored the actual generation 
of the turbulence, e.g., at the boundaries of the system, and the precise nature of the 
stirring force is not clear. In particular, they found that there was no a prion theory 
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for the exponents used to characterize the random-force autocorrelation function. 
First of all, for incompressible turbulence, 

~ + u(r, t) ■ Vu(r, t) - V[p/p\ + i/ 0 V 2 u(r, t) = F(r, t) (108) 

where the force F results from the interactions with the boundary and as such is 
not stochastic in nature. They took the view that a statistical force arises from the 
random force which represents the effective force felt at smaller length scales that 
results from turbulent, but deterministic, motion at larger scales as ‘transmitted’ by 
the nonlinear terms in NSE. They introduced a projection operator V such that 

u'(r, t) = ?u(r, t) (109) 

contains only high-wavenumber information. By applying V to Navier-Stokes equa- 
tion it follows 

~dt ~ + u '( r,t ^ ' Vu, ( r ’^ ~ ^[p'/p] + ^oV 2 u'(r ,t) = f(r, t) (110) 

where 

i(r,t) = VF (r,t) - P[u(r,<) • Vu(r,f)] + u'(r,t) • Vu'(r,f). (Ill) 

Yuan and Ronis (1992) noted that this ‘new’ force contains information about 
boundaries as well as the mode-coupling effects associated with velocity components 
on the injection scales. Since VF(r, t) = 0 away from boundaries, they defined the 
random stirring force used in RG studies to result from the mode coupling between the 
energy containing and the smaller scales. Since the motion on all scales is expected to 
be ‘chaotic’ including the energy containing scale, Yuan and Ronis (1992) expected 
that f(r, t) will have complicated, chaotic time and space dependences. Therefore, 
they identified this as the quantity which is actually modeled by a stochastic force in 
random stirring models of turbulence. 

The autocorrelation of the transverse parts of f(r, t) is assumed with a non-zero 
correlation time (colored noise). The major problem with such an assumption is that 
the resulting theory would not be invariant under Galilean transformation. Yuan and 
Ronis (1992) argued that there is no a priori reason why global Galilean invariance 
must hold. Their reason is that the random forces represents the effects of boundaries 
and these are not included in a Galilean transformation. Note that Eq. (110) is the 
subgrid NSE after the full fluctuating NSE is divided into the super- and sub- grid 
scales. However, the filtering operation is not used in the later development, and in 
particular, the path-integral RG works only with the full fluctuation Navier-Stokes 
equation. 
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7.4 The work of Carati 


The objective of Carati ( 1990a, b) is to modify the e-expansion of YO. 

Carati (1990a) first noticed the difference between the e- RG applications in Ma 
and Mazcnko (1975) and YO. In YO, the dimensionality d is considered as a fixed. 
A stochastic forcing is introduced to replace the initial and boundary conditions in 
Navier-Stokes equation. A correlation function of this forcing is assumed to following 
a powerlaw with exponent y. YO then defined y = e + d — 4 where e-expansion is 
performed for y. Carati (1990a) attempted to modify the YO procedure by introduc- 
ing a parametric dimension d = d 0 — 0 (e). But, as Carati (1990a) noted, the physical 
meaning of do is not clear a prior. 

The next effort of Carati (1990b) involved the extension of the white noise of 
the forcing correlation function to a colored one (again a colored correlation refers 
to a non-zero correlation in time). He argued, as did Yuan and Ronis (1992), that 
the correlation time as well as the correlation length can be important due to the 
appearance of macroscopic scale structures in turbulence. It is then possible to relate 
the expansion parameter to the stochastic forcing correlation. The result is that the 
parameter c is not treated as a small parameter and is fixed in its original physical 
value 4. The frequency correlation provides another free parameter through the pow- 
erlaw exponent q/3 where — 1 < 7 < 1 . The original YO’s e-expansion is replaced by 
arguing that the frequency integral is nearly divergent because of the assumed value 
of 7 . Alternatively, the e-expansion has now been replaced by a 7 -expansion where 
7=1 — 5. This ^-expansion was interpreted as a scheme in which a large amount of 
energy is injected into the system (Carati, 1990b). 

Although Carati viewed the colored forcing correlation as a good way to sidestep 
the e-expansion, he did not consider that the introduction of the colored noise violates 
the Galilean invariance. This important issue was pointed out by Yuan and Ronis 
(1992) in their analysis (sec previous subsection). 

7.5 The work of Liang and Diamond 

Liang and Diamond (1993) conducted a careful examination on the applicability of 
e-RG to two dimensional ( 2 D) flows. They examined both 2D fluid and 2 D Magnc- 
tohydrodynamicic (MHD) turbulence but our discussion will be limited to the fluid 
turbulence part of the paper only. 

The most important difference between 2 D and 3D turbulence is the direction of 
energy transfer. There is a dual-cascade phenomena, in which the energy is trans- 
fered from large to small scales, while enstrophy (defined as mean-square vorticity), 
is transferred from large scale to small scales (Kraichnan and Montgomery, 1980). 
Liang and Diamond (1993) found that when a proper inertial range energy spectrum 
is assumed, no fixed point can be found for positive or negative eddy viscosity. Therc- 
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fore, e-RG can not be used to analyze the physics of 2d turbulence. They found that 
the reason of this difficulty can be traced to the dual energy cascade. 

One important point should be stressed. Although Liang and Diamond (1993) 
followed most of the procedure of YO e-RG, they did not apply the e-expansion. 
Instead, they used the proper inertial range scaling as an input. 

To identify the reason behind the failure of e-RG in 2D, Liang and Diamond (1993) 
compared it to the EDQNM closure. In e-RG procedure, all calculations are one 
point and only eddy damping term appear. The nonlinear noise (nonlinear coupling) 
is replaced by adding stochastic forces. The major problem is that the properties of 
the stochastic forces can not be determined recursively. They are determined instead 
by the properties of the self-similarity range spectra, which in reality, are the goals 
of the calculation. On the other hand, a two-point closure calculation determines the 
form of the spectrum. Therefore, both the eddy damping and the nonlinear model 
coupling terms will now appear. Furthermore, while the e-RG only incorporate the 
nonlocal effects of the small scales on large scales, EDQNM (and DIA) take both 
effects into account, in principle. 

7.6 The work of Lam 

Lam’s approach (Lam 1992), although belonging to the FNS school of activity, is 
quite different from the other members of that group. It does not make use of the 
A-cxpansion procedure used by FNS, nor the correspondence principle and the e- 
expansion procedure used in the YO theory. His interpretation of e-RG is based on 
phenomenological approach. Zhou (1995) has shown that Lam’s model is essentially 
the physical space version of the classical closure theory (Leslie and Quarini, 1979) 
in spectral space and consider the corresponding treatment of the eddy viscosity and 
energy backseat ter. 

7.7 The work of Eyink 

Eyink (1994) used the momentum-shell RG method of Kadanoff- Wilson based on the 
Martin-Siggia-Rose (1978) field-theory formulation of stochastic dynamics. 

One of the major points in Eyink (1994) is that, contrary to the claim of YO, 
the higher-order nonlinear terms generated in their RG analysis are not irrelevant 
but marginal by power-counting. Because of this, the terms neglected by YO arc 
not necessarily small even for small e. The problem was traced to issues relevant to 
Galilean-invariant theories. It is shown not to occur in the FNS original analysis since 
they were dealing with fluctuation dynamics for equilibrium NS fluids. The issue of 
Galilean invariance has also been discussed by Zhou and Vahala (1993b). 

Eyink (1994) also presented his objections to the YO applications of e-RG to tur- 
bulence modeling. First, he repeated Kraichnan’s (1987) arguments that the power- 
laws derived by YO can be simply derived from dimensional analysis. Second, he 
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stated, as Zhou et al. (1988), that by setting e = 4, there is no reason to believe 
that additional nonlinearities arc negligible or insignificant to the Physics. Finally, 
the YO analysis does not include the rescaling procedure, which Eyink (1994) viewed 
as vital to error estimation. In fact, Eyink (1994) even stated that “YO thory is not 
an RG analysis at all!”. 
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Part IV 

RG-based turbulence modelling 

8 Modeling of Reynolds stress using RG 


Turbulent flows of practical interest involve a broad spectrum of length and time 
scales and require some type of modeling for Reynolds stresses. The most widely 
used arc the two-equation turbulence models based on the transport equations for 
parameters that involve the length and the time scales and require the simplest level 
of Reynolds stress closure that docs not depend specifically on the flow geometry. 
In its standard form the two-equation Reynolds stress turbulence models involve the 
turbulence kinetic energy and dissipation based on a Boussinesq type approximation 
involving isotropic eddy viscosity. Such a representation of turbulence is often not 
effective from both theoretical as well as phenomenological point of view. To overcome 
this, models that are nonlinear (i.e., quadratic) in the mean strain rate were proposed 
in the form of a constitutive relation (Speziale, 1987; Yoshizawa, 1984). These models 
arc capable of predicting the anisotropy in the Reynolds stresses in complex flows but 
require empirical evaluation of the model constants (Speziale, 1991). 

Two studies presented here 5 address the need for a more effective approach in the 
development of two-equation turbulence models, and in this context the renormaliza- 
tion group (RG) theory based models are considered (for a summary of the ICASE 
panel discussion on RG based turbulence modeling, see Zhou and Speziale (1994)). 
These models fall into two distinct categories: (a) e — RG, where a small parameter 
c is introduced into the exponent of the forcing correlation function (with the forcing 
function being introduced into the momentum equation), and the theory is then de- 
veloped for e -C 1 (Yakhot and Orszag, 1986); and (b) recursion- RG, which does not 
rely on an e-expansion, and treats explicitly the cubic nonlinearities introduced into 
the renormalized momentum equation (Zhou, et al. 1988). It should be noted that in 
the e- RG while all constants generated are evaluated in the limit t < 1, at the same 
time, all exponents that are e-dependent are evaluated at e = 4. In fact, e = 4 is re- 
quired in the e — RG to recover the Kolmogorov energy spectrum in the inertial range. 
In addition, e - RG theory can only take into account non-local interactions (Smith 
and Reynolds, 1992). On the other hand, recursion- /?G does not rely on an e-expan- 
sion, and treats explicitly the cubic nonlinearities introduced into the renormalized 
momentum equation. Moreover, recursion- i?C? can handle both local and non-local 
interactions. Effects such as the cusp behavior in the transport coefficients (Zhou 

5 Notc that recently Nagano and Itazu (1997) attempted to derived an eddy viscosity model using 
the iteractive averaging method 
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and Vahala, 1993a) which are due to local interactions and the cubic nonlinearities 
arc recovered in these theories. As in e — RG, the eddy viscosity is readily deter- 
mined from the solution of a relatively simple differential equation (Zhou and Vahala, 
1993a). However, unlike e — RG, the transport coefficients arc determined over the 
whole resolvable scales and not just in the wavenumber limit k — > 0. 

The contributions to the Reynolds stress tensor from the conventional time- 
averaging of the equations of motion are due to the nonlinear coupling term and its 
interaction with the velocity field (Zhou et al., 1993). In particular, the Reynolds 
stress tensor, t X] = + r~j~ , where the t£ + - part arises from the infrared limit 

of k — > 0 and is due to the uf — u1~ distant interaction limit while the - part 
arises from the 0 < k < k c spectrum and is due to the uf — uj local interaction 
limit. It should be noted that in e — RG approach, one takes the large-scale infrared 
limit k — > 0. In essence, this forces a spectral gap between the resolvable part of the 
flow field and the small unresolved scales. If this spectral gap were somehow present 
initially, it would be quickly populated in just a few eddy turn over times. Thus, 
retaining only the distant interactions may not be appropriate. In fact, it has been 
shown (Zhou and Vahala, 1993a) that the energy transfer function that corresponds to 
local interactions accounts for most of the energy flow out of the resolvable scales. It 
is, therefore, important to retain both local and nonlocal interactions in the modeling 
of the Reynolds stress and this can be readily achieved by recursion- RG. But in the 
e — RG model, the Reynolds stress t 13 = r+ + and is obtained purely from the uf — uj~ 
interaction in the small unresolved scale momentum equation. 

Rubinstein and Barton (1990) have derived a Reynolds stress model using the 
e — RG method (which corresponds to the infrared limit of k — > 0) of the following 
form: 

T+ + = + Vr(djUi + d t Uj) - ^[CrtidaUidaUjy 

+c r2 (d a u l d j u a + d a UjdiU a y + CrzidiUvdjU*)*]. ( 112 ) 

The constants C T \ = 0.034, C T 2 = 0.104 and C r3 = —0.014, and (...)* denotes 
the deviatoric part of the expression within the parenthesis. The first two terms 
correspond to the linear model and v T = C^K 2 je is the isotropic eddy- viscosity, 
where e is the turbulence dissipation and C ^ ~ 0.09 based on empirical data from 
equilibrium boundary layer flows. The above model is quadratic in mean strain rate, 
includes the effect of convection and diffusion and is qualitatively similar to other 
second order models (Speziale, 1991). 

Now, if one follows standard recursion- RG procedures (Zhou et ah, 1993), the 
relevant part of the small scale velocity field that contributes to is obtained 
in the wavenumber space. Transforming back to the physical space the following 
algebraic representation is obtained (Zhou et ah, 1993): 

T ij~( X ) - Cm^ldaUidpUidaUp + i~]} 
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(113) 


+C R ^[(d a U i )d j (d 1 Upd a d 0 U 1 ) +i~j] 

where £/* is the time-averaged mean velocity (and the second term inside the square 
brackets is obtained by switching the indices). However, for most flows of interest 
(inhomogeneous flows), the cut-off wavenumber k c varies considerably throughout the 
flow domain. Thus, if one wishes to use the local approximation, as well as retain 
k c ~ 0(1) effects so as to obtain an algebraic form for r^ < (x), then the coefficients 
Cr i and Cr 2 will be functions of the flow quantities. As a first attempt at applying this 
recursion- RG model we make the lowest order approximation that these coefficients 
are constants. In particular, for the present analysis, C R \ and C R 2 are taken to be 
0.025 and 0.342 x 10 -3 , respectively (Zhou et al., 1994). 

Combining (112) - (113) a formal expression for Reynolds stress which includes 
both the local and nonlocal interactions may be obtained. It is of some interest to 
note that integrity basis representations are commonly employed to represent the 
anisotropic part of the Reynolds stress tensor for three dimensional turbulent flows 
based on a systematic derivation from a hierarchy of second-order closure models 
(Gatski and Spczialc, 1993). It can readily be shown that the tensors that constitute 
the integrity basis are recovered for most part when the proposed recursion- RG model 
is recast appropriately (Zhou et al., 1994). 

The above expression for Reynolds stress (112) - (113) arc to be used along with 
the equations of motion by specifying turbulent kinetic energy and dissipation. In 
two-equation turbulence models, this closure is achieved through the development of 
transport equations for the turbulent kinetic energy and dissipation - quantities that 
are directly related to the length and time scales - of the following general form: 

d t K + UjdjK = V-£ + d t [(u 0 + Vr/a^K] (114) 

d t £ + UjdjE = CeiVe/K - C e2 e 2 jK + di[(u 0 + v T /a e )die\ (115) 

where, v — v G + u T is the total viscosity, V = ~r l j{dU i /dx ] ) is the turbulence 
production, e is the scalar turbulent dissipation rate. The quantities C e \, C e2 , 
are dimensionless and taken to be 1.44, 1.92, 1.0 and 1.3, respectively, consistent with 
the standard form of the two-equation K — e model (based on empirical data obtained 
from equilibrium boundary layer flows). We note that with various approximations, 
the e — RG-based formulations computed these constants as 1.42, 1.68, 0.719 and 
0.719, respectively (Yakhot and Orszag, 1986). 

The RG theory is utilized to develop Reynolds stress closure models for the pre- 
diction of turbulent separated flows. The combined model includes both the local and 
nonlocal interaction of all the relevant resolvable scales. The ability of the proposed 
model to accurately predict separated flows is analyzed from a combined theoretical 
and computational standpoint by considering turbulent flow past a backward facing 
step as a test case. The results obtained based on detailed computations demonstrate 
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that the RG model can yield very good predictions for the turbulent flow of an in- 
compressible viscous fluid over a backward-facing step (Zhou et al. , 1994). Thus, in 
spite of its well known deficiencies, when the anisotropy of the turbulent stresses arc 
properly accounted for, the two-equation turbulence models can be quite effective for 
the prediction of turbulent separated flows. 

9 RG based K — e model 

YO derived the K — e two equation model using e-RG method. Spczialc (1990) found 
that the original YO model performs poorly in homogeneous shear flow. The value C t \ 
derived in YO yields excessively large growth rate for the turbulent kinetic energy in 
homogeneous shear flow in comparision to both physical and numerical experiments 
(Spczialc, 1991). 

Smith and Reynolds (1992) found some algebraic error in the original derivation 
of YO. The coefficients of the dissipation term in e equation is not in good agreement 
with generally accepted values. Furthermore, YO’s derivation did not yield a term 
responsible for the production in the e equation. 

The original derivation of YO was revised by Yakhot and Smith (YS) (1992) by 
the following features: 

1. The ‘infrared cutoff’ of the random force, < rr >= 0 when 0 < k < A^ 

< fa(k, t)f 0 (k', t') >= D 0 k~ y D a p(k)8(k+k')6(t-t% A L < k < oo = 0. 0 < k < A L 

(116) 

This property is needed in the derivation of the equation for the mean rate of 
energy dissipation e (YS). 

2. The input of energy spectrum for the interval 0 < k < Al 

E{k) ~ k Q (117) 

is required to evaluate the integrals (with a = 2). 

3. Performing a Reynolds decomposition of T\ = — 2i^(VjUj)(Vjtq)(V/iq) into 
mean U and fluctuating u velocities. 

The derivation of YO and YS starts from dynamical equations for the homoge- 
neous part of the instantaneous rate of energy dissipation per unit mass e = ^o(Vjtq) 2 

Ve Ti 

Qs ' s / A v 

— +u l V i e = 2i/ 0 (VjWi)(Vj/j) — 2i/ 0 (VjUi)(VjU;)(ViUj) 
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T 2 

/ V 

- 2i/ 0 2 (V J V;u i ) 2 -2u 0 {V jUi){V iV jp) + v 0 VjVje (118) 

After some work, at the stirred fluids at the long-time and large-distance limit, 
the e-RG dissipation equation (YS, 1992) is found to be 


D t e = C t i(e/K)TijdjUi - C^e 2 /-^ + di(av T di£ ) -1Z (119) 


where C t \ = 1.42, C € 2 = 1.68 and 


TZ = 2u 0 Sij 


dui dui 
dxi dxj 


( 120 ) 


The e expansion procedure and above mentioned assumptions have been employed. 

The above e-RG dissipation equation is not closed. The neglect of TZ is a formally 
justified approximation at high Reynolds number if the hypothesis of local isotropy 
is invoked. Durbin and Speziale (1991) have questioned the validity of local isotropy 
in strongly strained turbulence flows. Yakhot et al. (1992) have proposed a model 
where TZ = TZ{rj), where the standard form of the model is recovered for TZ — * 0 in 
the limit of weak strains. Note that Durbin (1990) has already developed a model 
for the production of dissipation along these lines that was quadratic in the ratio 
of production to the dissipation and, hence, quartic in rj. Lam (1994) published a 
critique on the YS derivation. 

Iterating the expression for TZ using the Navier-Stokes equation will generate a 
power scries 


OO OIS 

1Z — u T S 3 jP r n ( )" (121) 

where S = (2S lJ S lJ ) l ' /2 . It is not possible to evalutc the summation since the values 
of coefficients arc unknown. 

The 1Z is modeled via three steps: 

1. The summation is performed for the geometric series for every three terms. 
This procedure reduces the numbers of unknown coefficients to one, j3. 

K° = v T s 3 ±(-gn—r = ( 122 ) 

2. Assuming that the fixed point value r/ 0 — 4.38 of the homogeneous shear flow 
in the equilibrium states, is invariant to dropping all terms but those in (122) above, 
Yakhot et al. (1992) postulated that 

v S 3 

^ = T+^ (1_77/?7o) - (123) 

3. One now further assumes that the isotropic Reynolds stress r VJ = —2 C^Kryj 
(% = SijK/s) 
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C^rfi 1 - r?/r? 0 ) e 3 rj(l - t//t/q) g 

1 + l3rf K 1 + /3rj 3 K ij ij 


(124) 


The undertermined constant (3 = 0.012 for the von Karman constant 0.4. The 
final e-RG dissipation rate transport equation is given 


D t e = C* a {ejK)rijdjUi - C e2 £ 2 /K + di(av T diE) 
where the coefficient C* A is given by 


(125) 


Q = c el 


vQ -v/vo) 

1 + /3rj 3 


(126) 


The model of Yakhot et al. (1992) has been tested for homogeneous shear flows 
and for flow over a backward facing step. Excellent results arc obtained in both cases. 
Recently, Yakhot and Orszag have extended the model and applied it to complex flows 
using the FLUENT code. 


Part V 

Conclusion 


In this review, we explained the concepts of terms renormalization and renormaliza- 
tion group by referencing to various physical systems, such as the ising model. We 
then present a comprehensive review on applications of the method of renormaliza- 
tion group to turbulence. These parts should be sufficient for readers who wishes to 
get a balanced view of RG in turbulence. For a few selected approaches, we have 
provided further technical details. We conclude with a discussion of the relevance 
and application of renormalization group to turbulence modelling. 
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